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

Merge branch 'state-integration-cleaning' into 'development' 

State integration cleaning

See merge request damask/DAMASK!155
This commit is contained in:
Franz Roters 2020-04-23 13:38:37 +02:00
parent 7452ebc727 84d6521183
commit 289808875c
18 changed files with 312 additions and 427 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
src/constitutive.f90
View File

@ -327,7 +327,7 @@ module constitutive
constitutive_initialFi, & constitutive_initialFi, &
constitutive_SandItsTangents, & constitutive_SandItsTangents, &
constitutive_collectDotState, & constitutive_collectDotState, &
constitutive_collectDeltaState, & constitutive_deltaState, &
constitutive_results constitutive_results
contains contains
@ -709,12 +709,14 @@ end subroutine constitutive_hooke_SandItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el, & !< element
phase, &
of
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
@ -730,16 +732,16 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip,
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
i, & !< counter in source loop i, & !< counter in source loop
instance, of instance
logical :: broken
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(phase)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
Mp = matmul(matmul(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) plasticityType: select case (phase_plasticity(phase))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState (Mp,instance,of) call plastic_isotropic_dotState (Mp,instance,of)
@ -760,10 +762,11 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip,
call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, &
instance,of,ip,el) instance,of,ip,el)
end select plasticityType end select plasticityType
broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of)))
SourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el)) SourceLoop: do i = 1, phase_Nsources(phase)
sourceType: select case (phase_source(i,material_phaseAt(ipc,el))) sourceType: select case (phase_source(i,phase))
case (SOURCE_damage_anisoBrittle_ID) sourceType case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress?
@ -775,25 +778,29 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip,
call source_damage_anisoDuctile_dotState ( ipc, ip, el) call source_damage_anisoDuctile_dotState ( ipc, ip, el)
case (SOURCE_thermal_externalheat_ID) sourceType case (SOURCE_thermal_externalheat_ID) sourceType
call source_thermal_externalheat_dotState(material_phaseAt(ipc,el),of) call source_thermal_externalheat_dotState(phase,of)
end select sourceType end select sourceType
broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of)))
enddo SourceLoop enddo SourceLoop
end subroutine constitutive_collectDotState end function constitutive_collectDotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> will return false if delta state is not needed/supported by the constitutive model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el, & !< element
phase, &
of
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola Kirchhoff stress S, & !< 2nd Piola Kirchhoff stress
Fe, & !< elastic deformation gradient Fe, & !< elastic deformation gradient
@ -802,35 +809,62 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
Mp Mp
integer :: & integer :: &
i, & i, &
instance, of instance, &
myOffset, &
mySize
logical :: &
broken
Mp = matmul(matmul(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(phase)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) plasticityType: select case (phase_plasticity(phase))
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_deltaState(Mp,instance,of) call plastic_kinehardening_deltaState(Mp,instance,of)
broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) call plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
case default
broken = .false.
end select plasticityType end select plasticityType
sourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el)) if(.not. broken) then
select case(phase_plasticity(phase))
case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
sourceType: select case (phase_source(i,material_phaseAt(ipc,el))) myOffset = plasticState(phase)%offsetDeltaState
mySize = plasticState(phase)%sizeDeltaState
plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = &
plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of)
end select
endif
sourceLoop: do i = 1, phase_Nsources(phase)
sourceType: select case (phase_source(i,phase))
case (SOURCE_damage_isoBrittle_ID) sourceType case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, &
ipc, ip, el) ipc, ip, el)
broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of)))
if(.not. broken) then
myOffset = sourceState(phase)%p(i)%offsetDeltaState
mySize = sourceState(phase)%p(i)%sizeDeltaState
sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) = &
sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) + sourceState(phase)%p(i)%deltaState(1:mySize,of)
endif
end select sourceType end select sourceType
enddo SourceLoop enddo SourceLoop
end subroutine constitutive_collectDeltaState end function constitutive_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

2
src/constitutive_plastic_disloUCLA.f90
View File

@ -209,7 +209,7 @@ module subroutine plastic_disloUCLA_init
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

2
src/constitutive_plastic_dislotwin.f90
View File

@ -399,7 +399,7 @@ module subroutine plastic_dislotwin_init
+ size(['f_tr']) * prm%sum_N_tr + size(['f_tr']) * prm%sum_N_tr
sizeState = sizeDotState sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! locally defined state aliases and initialization of state0 and atol

2
src/constitutive_plastic_isotropic.f90
View File

@ -117,7 +117,7 @@ module subroutine plastic_isotropic_init
sizeDotState = size(['xi ','accumulated_shear']) sizeDotState = size(['xi ','accumulated_shear'])
sizeState = sizeDotState sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

2
src/constitutive_plastic_kinehardening.f90
View File

@ -164,7 +164,7 @@ module subroutine plastic_kinehardening_init
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

2
src/constitutive_plastic_none.f90
View File

@ -29,7 +29,7 @@ module subroutine plastic_none_init
if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
call material_allocatePlasticState(p,NipcMyPhase,0,0,0) call material_allocateState(plasticState(p),NipcMyPhase,0,0,0)
enddo enddo

5
src/constitutive_plastic_nonlocal.f90
View File

@ -320,6 +320,7 @@ module subroutine plastic_nonlocal_init
prm%fEdgeMultiplication = config%getFloat('edgemultiplication') prm%fEdgeMultiplication = config%getFloat('edgemultiplication')
prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/') prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity checks
if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' burgers' if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' burgers'
@ -384,9 +385,9 @@ module subroutine plastic_nonlocal_init
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState sizeDeltaState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
plasticState(p)%nonlocal = .true. plasticState(p)%nonlocal = config%KeyExists('/nonlocal/')
plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)

2
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -213,7 +213,7 @@ module subroutine plastic_phenopowerlaw_init
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

544
src/crystallite.f90
View File

@ -15,7 +15,6 @@ module crystallite
use DAMASK_interface use DAMASK_interface
use config use config
use debug use debug
use numerics
use rotations use rotations
use math use math
use FEsolving use FEsolving
@ -70,9 +69,7 @@ module crystallite
logical, dimension(:,:,:), allocatable, public :: & logical, dimension(:,:,:), allocatable, public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
logical, dimension(:,:,:), allocatable :: & logical, dimension(:,:,:), allocatable :: &
crystallite_converged, & !< convergence flag crystallite_converged !< convergence flag
crystallite_todo, & !< flag to indicate need for further computation
crystallite_localPlasticity !< indicates this grain to have purely local constitutive law
type :: tOutput !< new requested output (per phase) type :: tOutput !< new requested output (per phase)
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
@ -84,7 +81,8 @@ module crystallite
integer :: & integer :: &
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
nState, & !< state loop limit nState, & !< state loop limit
nStress !< stress loop limit nStress, & !< stress loop limit
integrator !< integration scheme (ToDo: better use a string)
real(pReal) :: & real(pReal) :: &
subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback
subStepSizeCryst, & !< size of first substep when cutback subStepSizeCryst, & !< size of first substep when cutback
@ -98,7 +96,7 @@ module crystallite
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
procedure(), pointer :: integrateState procedure(integrateStateFPI), pointer :: integrateState
public :: & public :: &
crystallite_init, & crystallite_init, &
@ -159,9 +157,7 @@ subroutine crystallite_init
allocate(crystallite_orientation(cMax,iMax,eMax)) allocate(crystallite_orientation(cMax,iMax,eMax))
allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.)
allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.)
allocate(crystallite_todo(cMax,iMax,eMax), source=.false.)
allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) allocate(crystallite_converged(cMax,iMax,eMax), source=.true.)
num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal) num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal)
@ -177,6 +173,8 @@ subroutine crystallite_init
num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1) num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1)
num%integrator = config_numerics%getInt ('integrator', defaultVal=1)
num%nState = config_numerics%getInt ('nstate', defaultVal=20) num%nState = config_numerics%getInt ('nstate', defaultVal=20)
num%nStress = config_numerics%getInt ('nstress', defaultVal=40) num%nStress = config_numerics%getInt ('nstress', defaultVal=40)
@ -193,10 +191,14 @@ subroutine crystallite_init
if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum')
if(num%integrator < 1 .or. num%integrator > 5) &
call IO_error(301,ext_msg='integrator')
if(num%nState < 1) call IO_error(301,ext_msg='nState') if(num%nState < 1) call IO_error(301,ext_msg='nState')
if(num%nStress< 1) call IO_error(301,ext_msg='nStress') if(num%nStress< 1) call IO_error(301,ext_msg='nStress')
select case(numerics_integrator)
select case(num%integrator)
case(1) case(1)
integrateState => integrateStateFPI integrateState => integrateStateFPI
case(2) case(2)
@ -234,7 +236,6 @@ subroutine crystallite_init
/ math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal)
crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e) crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e)
crystallite_F0(1:3,1:3,c,i,e) = math_I3 crystallite_F0(1:3,1:3,c,i,e) = math_I3
crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phaseAt(c,e))
crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), & crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), &
crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration
crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e)
@ -244,7 +245,7 @@ subroutine crystallite_init
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal?
crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedFi0 = crystallite_Fi0 crystallite_partionedFi0 = crystallite_Fi0
@ -271,9 +272,8 @@ subroutine crystallite_init
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then
write(6,'(a42,1x,i10)') ' # of elements: ', eMax write(6,'(a42,1x,i10)') ' # of elements: ', eMax
write(6,'(a42,1x,i10)') 'max # of integration points/element: ', iMax write(6,'(a42,1x,i10)') ' # of integration points/element: ', iMax
write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax
write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity)
flush(6) flush(6)
endif endif
@ -301,6 +301,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
e, & !< counter in element loop e, & !< counter in element loop
startIP, endIP, & startIP, endIP, &
s s
logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 &
@ -344,7 +345,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e)
crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subFrac(c,i,e) = 0.0_pReal
crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst
crystallite_todo(c,i,e) = .true. todo(c,i,e) = .true.
crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst
endif homogenizationRequestsCalculation endif homogenizationRequestsCalculation
enddo; enddo enddo; enddo
@ -361,7 +362,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
endif singleRun endif singleRun
NiterationCrystallite = 0 NiterationCrystallite = 0
cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) cutbackLooping: do while (any(todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2))))
NiterationCrystallite = NiterationCrystallite + 1 NiterationCrystallite = NiterationCrystallite + 1
#ifdef DEBUG #ifdef DEBUG
@ -380,8 +381,8 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), &
num%stepIncreaseCryst * crystallite_subStep(c,i,e)) num%stepIncreaseCryst * crystallite_subStep(c,i,e))
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on?
if (crystallite_todo(c,i,e)) then if (todo(c,i,e)) then
crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e)
crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e)
crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e)
@ -415,12 +416,12 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
enddo enddo
! cant restore dotState here, since not yet calculated in first cutback after initialization ! cant restore dotState here, since not yet calculated in first cutback after initialization
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! prepare for integration ! prepare for integration
if (crystallite_todo(c,i,e)) then if (todo(c,i,e)) then
crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) &
+ crystallite_subStep(c,i,e) *( crystallite_partionedF (1:3,1:3,c,i,e) & + crystallite_subStep(c,i,e) *( crystallite_partionedF (1:3,1:3,c,i,e) &
-crystallite_partionedF0(1:3,1:3,c,i,e)) -crystallite_partionedF0(1:3,1:3,c,i,e))
@ -438,9 +439,9 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! integrate --- requires fully defined state array (basic + dependent state) ! integrate --- requires fully defined state array (basic + dependent state)
if (any(crystallite_todo)) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation if (any(todo)) call integrateState(todo) ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation
where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further
crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
enddo cutbackLooping enddo cutbackLooping
@ -610,14 +611,16 @@ subroutine crystallite_orientations
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
nonlocalPresent: if (any(plasticState%nonLocal)) then nonlocalPresent: if (any(plasticState%nonlocal)) then
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) if (plasticState(material_phaseAt(1,e))%nonlocal) then
if (plasticState(material_phaseAt(1,e))%nonLocal) & do i = FEsolving_execIP(1),FEsolving_execIP(2)
call plastic_nonlocal_updateCompatibility(crystallite_orientation, & call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(i,e)),i,e) phase_plasticityInstance(material_phaseAt(i,e)),i,e)
enddo; enddo enddo
endif
enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
endif nonlocalPresent endif nonlocalPresent
@ -777,7 +780,7 @@ end subroutine crystallite_results
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction !> intermediate acceleration of the Newton-Raphson correction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical function integrateStress(ipc,ip,el,timeFraction) function integrateStress(ipc,ip,el,timeFraction) result(broken)
integer, intent(in):: el, & ! element index integer, intent(in):: el, & ! element index
ip, & ! integration point index ip, & ! integration point index
@ -834,9 +837,9 @@ logical function integrateStress(ipc,ip,el,timeFraction)
p, & p, &
jacoCounterLp, & jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update jacoCounterLi ! counters to check for Jacobian update
logical :: error logical :: error,broken
integrateStress = .false. broken = .true.
if (present(timeFraction)) then if (present(timeFraction)) then
dt = crystallite_subdt(ipc,ip,el) * timeFraction dt = crystallite_subdt(ipc,ip,el) * timeFraction
@ -847,6 +850,9 @@ logical function integrateStress(ipc,ip,el,timeFraction)
F = crystallite_subF(1:3,1:3,ipc,ip,el) F = crystallite_subF(1:3,1:3,ipc,ip,el)
endif endif
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,ipc,ip,el), &
crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el)
Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess
Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess
@ -977,7 +983,6 @@ logical function integrateStress(ipc,ip,el,timeFraction)
call math_invert33(Fp_new,devNull,error,invFp_new) call math_invert33(Fp_new,devNull,error,invFp_new)
if (error) return ! error if (error) return ! error
integrateStress = .true.
crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new)))
crystallite_S (1:3,1:3,ipc,ip,el) = S crystallite_S (1:3,1:3,ipc,ip,el) = S
crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess
@ -985,6 +990,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize
crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new
crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new)
broken = .false.
end function integrateStress end function integrateStress
@ -993,8 +999,9 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method !> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize !> using Fixed Point Iteration to adapt the stepsize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateFPI subroutine integrateStateFPI(todo)
logical, dimension(:,:,:), intent(in) :: todo
integer :: & integer :: &
NiterationState, & !< number of iterations in state loop NiterationState, & !< number of iterations in state loop
e, & !< element index in element loop e, & !< element index in element loop
@ -1003,118 +1010,107 @@ subroutine integrateStateFPI
p, & p, &
c, & c, &
s, & s, &
sizeDotState size_pl
integer, dimension(maxval(phase_Nsources)) :: &
size_so
real(pReal) :: & real(pReal) :: &
zeta zeta
real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: &
r ! state residuum r ! state residuum
real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: &
plastic_dotState
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState) !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState,source_dotState,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
do s = 1, phase_Nsources(p) if(broken) cycle
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState size_pl = plasticState(p)%sizeDotState
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) &
+ plasticState(p)%dotState (1:sizeDotState,c) & + plasticState(p)%dotState (1:size_pl,c) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
plastic_dotState_p2 = 0.0_pReal * plasticState(p)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer plastic_dotState(1:size_pl,2) = 0.0_pReal
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState size_so(s) = sourceState(p)%p(s)%sizeDotState
sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) &
+ sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:size_so(s),c) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
source_dotState(1:sizeDotState,2,s) = 0.0_pReal source_dotState(1:size_so(s),2,s) = 0.0_pReal
enddo enddo
iteration: do NiterationState = 1, num%nState iteration: do NiterationState = 1, num%nState
if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1)
plastic_dotState_p1 = plasticState(p)%dotState(:,c) plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c)
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s)
if(nIterationState > 1) source_dotState(1:sizeDotState,2,s) = source_dotState(1:sizeDotState,1,s) source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c)
source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c)
enddo enddo
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e)
crystallite_Fp(1:3,1:3,g,i,e), & if(broken) exit iteration
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. crystallite_todo(g,i,e)) exit iteration crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) exit iteration
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),&
crystallite_partionedF0, & plastic_dotState(1:size_pl,2))
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. crystallite_todo(g,i,e)) exit iteration
sizeDotState = plasticState(p)%sizeDotState
zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2)
plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta &
+ plastic_dotState_p1 * (1.0_pReal - zeta) + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta)
r(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & r(1:size_pl) = plasticState(p)%state (1:size_pl,c) &
- plasticState(p)%subState0(1:sizeDotState,c) & - plasticState(p)%subState0(1:size_pl,c) &
- plasticState(p)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e)
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) &
- r(1:sizeDotState) - r(1:size_pl)
crystallite_converged(g,i,e) = converged(r(1:sizeDotState), & crystallite_converged(g,i,e) = converged(r(1:size_pl), &
plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%state(1:size_pl,c), &
plasticState(p)%atol(1:sizeDotState)) plasticState(p)%atol(1:size_pl))
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
zeta = damper(sourceState(p)%p(s)%dotState(:,c), & zeta = damper(sourceState(p)%p(s)%dotState(:,c), &
source_dotState(1:sizeDotState,1,s),& source_dotState(1:size_so(s),1,s),&
source_dotState(1:sizeDotState,2,s)) source_dotState(1:size_so(s),2,s))
sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta &
+ source_dotState(1:sizeDotState,1,s)* (1.0_pReal - zeta) + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta)
r(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) & r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) &
- sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - sourceState(p)%p(s)%subState0(1:size_so(s),c) &
- sourceState(p)%p(s)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e)
sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) & sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) &
- r(1:sizeDotState) - r(1:size_so(s))
crystallite_converged(g,i,e) = & crystallite_converged(g,i,e) = &
crystallite_converged(g,i,e) .and. converged(r(1:sizeDotState), & crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), &
sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%state(1:size_so(s),c), &
sourceState(p)%p(s)%atol(1:sizeDotState)) sourceState(p)%p(s)%atol(1:size_so(s)))
enddo enddo
if(crystallite_converged(g,i,e)) then if(crystallite_converged(g,i,e)) then
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
exit iteration exit iteration
endif endif
enddo iteration enddo iteration
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
nonlocalBroken = .true.
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -1149,7 +1145,9 @@ end subroutine integrateStateFPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method !> @brief integrate state with 1st order explicit Euler method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateEuler subroutine integrateStateEuler(todo)
logical, dimension(:,:,:), intent(in) :: todo
integer :: & integer :: &
e, & !< element index in element loop e, & !< element index in element loop
@ -1160,29 +1158,25 @@ subroutine integrateStateEuler
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE (sizeDotState,p,c) !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
do s = 1, phase_Nsources(p) if(broken) cycle
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
@ -1195,21 +1189,15 @@ subroutine integrateStateEuler
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
enddo enddo
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & crystallite_Fe(1:3,1:3,g,i,e), &
nonlocalBroken = .true. crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(.not. crystallite_todo(g,i,e)) cycle if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
if(broken) cycle
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
crystallite_converged(g,i,e) = crystallite_todo(g,i,e)
broken = integrateStress(g,i,e)
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
crystallite_converged(g,i,e) = .not. broken
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -1222,7 +1210,9 @@ end subroutine integrateStateEuler
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size !> @brief integrate stress, state with 1st order Euler method with adaptive step size
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateAdaptiveEuler subroutine integrateStateAdaptiveEuler(todo)
logical, dimension(:,:,:), intent(in) :: todo
integer :: & integer :: &
e, & ! element index in element loop e, & ! element index in element loop
@ -1233,32 +1223,28 @@ subroutine integrateStateAdaptiveEuler
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic
real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then broken = .false.
p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) cycle
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1274,36 +1260,23 @@ subroutine integrateStateAdaptiveEuler
+ sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e)
enddo enddo
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & crystallite_Fe(1:3,1:3,g,i,e), &
nonlocalBroken = .true. crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(.not. crystallite_todo(g,i,e)) cycle if(broken) cycle
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e)
crystallite_Fp(1:3,1:3,g,i,e), & if(broken) cycle
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) cycle
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) &
+ 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), &
plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%state(1:sizeDotState,c), &
@ -1311,7 +1284,6 @@ subroutine integrateStateAdaptiveEuler
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState sizeDotState = sourceState(p)%p(s)%sizeDotState
crystallite_converged(g,i,e) = & crystallite_converged(g,i,e) = &
crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) &
+ 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), &
@ -1320,6 +1292,7 @@ subroutine integrateStateAdaptiveEuler
enddo enddo
endif endif
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -1331,18 +1304,20 @@ end subroutine integrateStateAdaptiveEuler
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 4th order explicit Runge Kutta method !> @brief integrate stress, state with 4th order explicit Runge Kutta method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateRK4 subroutine integrateStateRK4(todo)
real(pReal), dimension(3,3), parameter :: & logical, dimension(:,:,:), intent(in) :: todo
A = reshape([&
real(pReal), dimension(3,3), parameter :: &
A = reshape([&
0.5_pReal, 0.0_pReal, 0.0_pReal, & 0.5_pReal, 0.0_pReal, 0.0_pReal, &
0.0_pReal, 0.5_pReal, 0.0_pReal, & 0.0_pReal, 0.5_pReal, 0.0_pReal, &
0.0_pReal, 0.0_pReal, 1.0_pReal], & 0.0_pReal, 0.0_pReal, 1.0_pReal], &
[3,3]) [3,3])
real(pReal), dimension(3), parameter :: & real(pReal), dimension(3), parameter :: &
CC = [0.5_pReal, 0.5_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration CC = [0.5_pReal, 0.5_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
real(pReal), dimension(4), parameter :: & real(pReal), dimension(4), parameter :: &
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] ! weight of slope used for Runge Kutta integration (final weight divided by 6) B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] ! weight of slope used for Runge Kutta integration (final weight divided by 6)
integer :: & integer :: &
e, & ! element index in element loop e, & ! element index in element loop
@ -1355,31 +1330,28 @@ subroutine integrateStateRK4
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState
real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RK4dotState,source_RK4dotState) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then broken = .false.
p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) cycle
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
do stage = 1,3 do stage = 1,3
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1413,31 +1385,18 @@ subroutine integrateStateRK4
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
enddo enddo
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e,CC(stage))
crystallite_Fp(1:3,1:3,g,i,e), & if(broken) exit
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e)*CC(stage), g,i,e) crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) exit
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
enddo enddo
if(broken) cycle
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1459,25 +1418,16 @@ subroutine integrateStateRK4
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
enddo enddo
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & crystallite_Fe(1:3,1:3,g,i,e), &
nonlocalBroken = .true. crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(.not. crystallite_todo(g,i,e)) cycle if(broken) cycle
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e)
crystallite_Fp(1:3,1:3,g,i,e), & crystallite_converged(g,i,e) = .not. broken
g, i, e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
crystallite_todo(g,i,e) = integrateStress(g,i,e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken
endif endif
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -1490,7 +1440,9 @@ end subroutine integrateStateRK4
!> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with !> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with
!> adaptive step size (use 5th order solution to advance = "local extrapolation") !> adaptive step size (use 5th order solution to advance = "local extrapolation")
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45 subroutine integrateStateRKCK45(todo)
logical, dimension(:,:,:), intent(in) :: todo
real(pReal), dimension(5,5), parameter :: & real(pReal), dimension(5,5), parameter :: &
A = reshape([& A = reshape([&
@ -1523,31 +1475,27 @@ subroutine integrateStateRKCK45
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState
real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then broken = .false.
p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) cycle
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
do stage = 1,5 do stage = 1,5
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1581,31 +1529,18 @@ subroutine integrateStateRKCK45
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
enddo enddo
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e,CC(stage))
crystallite_Fp(1:3,1:3,g,i,e), & if(broken) exit
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e)*CC(stage), g,i,e) crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) exit
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
enddo enddo
if(broken) cycle
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1614,7 +1549,7 @@ subroutine integrateStateRKCK45
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+ plasticState(p)%dotState (1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
crystallite_todo(g,i,e) = converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) & broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) &
* crystallite_subdt(g,i,e), & * crystallite_subdt(g,i,e), &
plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%atol(1:sizeDotState)) plasticState(p)%atol(1:sizeDotState))
@ -1627,35 +1562,28 @@ subroutine integrateStateRKCK45
sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) &
+ sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & broken = broken .and. .not. &
converged(matmul(source_RKdotState(1:sizeDotState,1:6,s),DB) & converged(matmul(source_RKdotState(1:sizeDotState,1:6,s),DB) &
* crystallite_subdt(g,i,e), & * crystallite_subdt(g,i,e), &
sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%state(1:sizeDotState,c), &
sourceState(p)%p(s)%atol(1:sizeDotState)) sourceState(p)%p(s)%atol(1:sizeDotState))
enddo enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & if(broken) cycle
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & crystallite_Fe(1:3,1:3,g,i,e), &
nonlocalBroken = .true. crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(.not. crystallite_todo(g,i,e)) cycle if(broken) cycle
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e)
crystallite_Fp(1:3,1:3,g,i,e), & crystallite_converged(g,i,e) = .not. broken
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken
endif endif
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
if (nonlocalBroken) call nonlocalConvergenceCheck if(nonlocalBroken) call nonlocalConvergenceCheck
end subroutine integrateStateRKCK45 end subroutine integrateStateRKCK45
@ -1666,7 +1594,16 @@ end subroutine integrateStateRKCK45
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine nonlocalConvergenceCheck subroutine nonlocalConvergenceCheck
where( .not. crystallite_localPlasticity) crystallite_converged = .false. integer :: e,i,p
!$OMP PARALLEL DO PRIVATE(p)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
p = material_phaseAt(1,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
if(plasticState(p)%nonlocal) crystallite_converged(1,i,e) = .false.
enddo
enddo
!$OMP END PARALLEL DO
end subroutine nonlocalConvergenceCheck end subroutine nonlocalConvergenceCheck
@ -1688,59 +1625,6 @@ logical pure function converged(residuum,state,atol)
end function converged end function converged
!--------------------------------------------------------------------------------------------------
!> @brief calculates a jump in the state according to the current state and the current stress
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!--------------------------------------------------------------------------------------------------
logical function stateJump(ipc,ip,el)
integer, intent(in):: &
el, & ! element index
ip, & ! integration point index
ipc ! grain index
integer :: &
c, &
p, &
mySource, &
myOffset, &
mySize
c = material_phaseMemberAt(ipc,ip,el)
p = material_phaseAt(ipc,el)
call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), &
crystallite_Fe(1:3,1:3,ipc,ip,el), &
crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc,ip,el)
myOffset = plasticState(p)%offsetDeltaState
mySize = plasticState(p)%sizeDeltaState
if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then
stateJump = .false.
return
endif
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = &
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c)
do mySource = 1, phase_Nsources(p)
myOffset = sourceState(p)%p(mySource)%offsetDeltaState
mySize = sourceState(p)%p(mySource)%sizeDeltaState
if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then
stateJump = .false.
return
endif
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = &
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c)
enddo
stateJump = .true.
end function stateJump
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current restart information (Field and constitutive data) to file. !> @brief Write current restart information (Field and constitutive data) to file.
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively

66
src/material.f90
View File

@ -11,7 +11,6 @@ module material
use results use results
use IO use IO
use debug use debug
use numerics
use rotations use rotations
use discretization use discretization
@ -174,8 +173,7 @@ module material
public :: & public :: &
material_init, & material_init, &
material_allocatePlasticState, & material_allocateState, &
material_allocateSourceState, &
ELASTICITY_HOOKE_ID ,& ELASTICITY_HOOKE_ID ,&
PLASTICITY_NONE_ID, & PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, & PLASTICITY_ISOTROPIC_ID, &
@ -700,63 +698,35 @@ end subroutine material_parseTexture
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates the plastic state of a phase !> @brief Allocate the components of the state structure for a given phase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine material_allocatePlasticState(phase,NipcMyPhase,& subroutine material_allocateState(state, &
sizeState,sizeDotState,sizeDeltaState) NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
class(tState), intent(out) :: &
state
integer, intent(in) :: & integer, intent(in) :: &
phase, &
NipcMyPhase, & NipcMyPhase, &
sizeState, & sizeState, &
sizeDotState, & sizeDotState, &
sizeDeltaState sizeDeltaState
plasticState(phase)%sizeState = sizeState state%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState state%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState state%sizeDeltaState = sizeDeltaState
plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(plasticState(phase)%atol (sizeState), source=0.0_pReal) allocate(state%atol (sizeState), source=0.0_pReal)
allocate(plasticState(phase)%state0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%state0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%partionedState0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%partionedState0(sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%subState0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%state (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(state%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) allocate(state%deltaState(sizeDeltaState,NipcMyPhase), source=0.0_pReal)
end subroutine material_allocatePlasticState end subroutine material_allocateState
!--------------------------------------------------------------------------------------------------
!> @brief allocates the source state of a phase
!--------------------------------------------------------------------------------------------------
subroutine material_allocateSourceState(phase,of,NipcMyPhase,&
sizeState,sizeDotState,sizeDeltaState)
integer, intent(in) :: &
phase, &
of, &
NipcMyPhase, &
sizeState, sizeDotState,sizeDeltaState
sourceState(phase)%p(of)%sizeState = sizeState
sourceState(phase)%p(of)%sizeDotState = sizeDotState
sourceState(phase)%p(of)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(of)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(sourceState(phase)%p(of)%atol (sizeState), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%state0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%partionedState0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%subState0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
end subroutine material_allocateSourceState
end module material end module material

23
src/numerics.f90
View File

@ -20,8 +20,7 @@ module numerics
iJacoStiffness = 1, & !< frequency of stiffness update iJacoStiffness = 1, & !< frequency of stiffness update
randomSeed = 0, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed randomSeed = 0, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed
worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only)
worldsize = 1, & !< MPI worldsize (/=1 for MPI simulations only) worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only)
numerics_integrator = 1 !< method used for state integration Default 1: fix-point iteration
integer(4), protected, public :: & integer(4), protected, public :: &
DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
real(pReal), protected, public :: & real(pReal), protected, public :: &
@ -134,8 +133,6 @@ subroutine numerics_init
defgradTolerance = IO_floatValue(line,chunkPos,2) defgradTolerance = IO_floatValue(line,chunkPos,2)
case ('ijacostiffness') case ('ijacostiffness')
iJacoStiffness = IO_intValue(line,chunkPos,2) iJacoStiffness = IO_intValue(line,chunkPos,2)
case ('integrator')
numerics_integrator = IO_intValue(line,chunkPos,2)
case ('usepingpong') case ('usepingpong')
usepingpong = IO_intValue(line,chunkPos,2) > 0 usepingpong = IO_intValue(line,chunkPos,2) > 0
case ('unitlength') case ('unitlength')
@ -176,6 +173,11 @@ subroutine numerics_init
case ('maxstaggerediter') case ('maxstaggerediter')
stagItMax = IO_intValue(line,chunkPos,2) stagItMax = IO_intValue(line,chunkPos,2)
#ifdef PETSC
case ('petsc_options')
petsc_options = trim(line(chunkPos(4):))
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! spectral parameters ! spectral parameters
#ifdef Grid #ifdef Grid
@ -187,8 +189,6 @@ subroutine numerics_init
err_stress_tolrel = IO_floatValue(line,chunkPos,2) err_stress_tolrel = IO_floatValue(line,chunkPos,2)
case ('err_stress_tolabs') case ('err_stress_tolabs')
err_stress_tolabs = IO_floatValue(line,chunkPos,2) err_stress_tolabs = IO_floatValue(line,chunkPos,2)
case ('petsc_options')
petsc_options = trim(line(chunkPos(4):))
case ('err_curl_tolabs') case ('err_curl_tolabs')
err_curl_tolAbs = IO_floatValue(line,chunkPos,2) err_curl_tolAbs = IO_floatValue(line,chunkPos,2)
case ('err_curl_tolrel') case ('err_curl_tolrel')
@ -206,8 +206,6 @@ subroutine numerics_init
integrationorder = IO_intValue(line,chunkPos,2) integrationorder = IO_intValue(line,chunkPos,2)
case ('structorder') case ('structorder')
structorder = IO_intValue(line,chunkPos,2) structorder = IO_intValue(line,chunkPos,2)
case ('petsc_options')
petsc_options = trim(line(chunkPos(4):))
case ('bbarstabilisation') case ('bbarstabilisation')
BBarStabilisation = IO_intValue(line,chunkPos,2) > 0 BBarStabilisation = IO_intValue(line,chunkPos,2) > 0
#endif #endif
@ -223,7 +221,6 @@ subroutine numerics_init
! writing parameters to output ! writing parameters to output
write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,1x,i8)') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength
@ -266,7 +263,6 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel
write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha
write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options)
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -274,16 +270,17 @@ subroutine numerics_init
#ifdef FEM #ifdef FEM
write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder
write(6,'(a24,1x,i8)') ' structOrder: ',structOrder write(6,'(a24,1x,i8)') ' structOrder: ',structOrder
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options)
write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation
#endif #endif
#ifdef PETSC
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options)
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity checks
if (defgradTolerance <= 0.0_pReal) call IO_error(301,ext_msg='defgradTolerance') if (defgradTolerance <= 0.0_pReal) call IO_error(301,ext_msg='defgradTolerance')
if (iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness') if (iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness')
if (numerics_integrator <= 0 .or. numerics_integrator >= 6) &
call IO_error(301,ext_msg='integrator')
if (numerics_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength') if (numerics_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength')
if (residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
if (itmax <= 1) call IO_error(301,ext_msg='itmax') if (itmax <= 1) call IO_error(301,ext_msg='itmax')

3
src/prec.f90
View File

@ -53,8 +53,7 @@ module prec
logical :: & logical :: &
nonlocal = .false. nonlocal = .false.
real(pReal), pointer, dimension(:,:) :: & real(pReal), pointer, dimension(:,:) :: &
slipRate, & !< slip rate slipRate !< slip rate
accumulatedSlip !< accumulated plastic slip
end type end type
type :: tSourceState type :: tSourceState

2
src/source_damage_anisoBrittle.f90
View File

@ -107,7 +107,7 @@ subroutine source_damage_anisoBrittle_init
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critDisp' if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critDisp'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'

2
src/source_damage_anisoDuctile.f90
View File

@ -89,7 +89,7 @@ subroutine source_damage_anisoDuctile_init
if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain' if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisoductile_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisoductile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'

2
src/source_damage_isoBrittle.f90
View File

@ -83,7 +83,7 @@ subroutine source_damage_isoBrittle_init
if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy' if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,1) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,1)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('isobrittle_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = config%getFloat('isobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'

2
src/source_damage_isoDuctile.f90
View File

@ -82,7 +82,7 @@ subroutine source_damage_isoDuctile_init
if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain' if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('isoductile_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = config%getFloat('isoductile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'

2
src/source_thermal_dissipation.f90
View File

@ -67,7 +67,7 @@ subroutine source_thermal_dissipation_init
prm%kappa = config%getFloat('dissipation_coldworkcoeff') prm%kappa = config%getFloat('dissipation_coldworkcoeff')
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,0,0,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,0,0,0)
end associate end associate
enddo enddo

2
src/source_thermal_externalheat.f90
View File

@ -74,7 +74,7 @@ subroutine source_thermal_externalheat_init
prm%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(prm%time)) prm%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(prm%time))
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
end associate end associate
enddo enddo