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

Merge branch 'development' into DADF5-multiprocessing

This commit is contained in:
Martin Diehl 2020-02-25 17:58:07 +01:00
parent 8e6ce6cbd6 822e6b7199
commit 75528064e9
13 changed files with 740 additions and 1137 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
v2.0.3-1732-gab88ffd2
v2.0.3-1747-ga2e8e5b1

9
python/damask/dadf5.py
View File

@ -111,10 +111,13 @@ class DADF5():
if action == 'set':
self.selection[what] = valid
elif action == 'add':
self.selection[what] = list(existing.union(valid))
add=existing.union(valid)
add_sorted=sorted(add, key=lambda x: int("".join([i for i in x if i.isdigit()])))
self.selection[what] = add_sorted
elif action == 'del':
self.selection[what] = list(existing.difference_update(valid))
diff=existing.difference(valid)
diff_sorted=sorted(diff, key=lambda x: int("".join([i for i in x if i.isdigit()])))
self.selection[what] = diff_sorted
def __time_to_inc(self,start,end):
selected = []

255
src/constitutive_plastic_disloUCLA.f90
View File

@ -6,21 +6,10 @@
!> @brief crystal plasticity model for bcc metals, especially Tungsten
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_disloUCLA
real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
enum, bind(c)
enumerator :: &
undefined_ID, &
rho_mob_ID, &
rho_dip_ID, &
dot_gamma_sl_ID, &
gamma_sl_ID, &
Lambda_sl_ID, &
tau_pass_ID
end enum
type :: tParameters
real(pReal) :: &
aTol_rho, &
@ -28,7 +17,7 @@ submodule(constitutive) plastic_disloUCLA
mu, &
D_0, & !< prefactor for self-diffusion coefficient
Q_cl !< activation energy for dislocation climb
real(pReal), dimension(:), allocatable :: &
real(pReal), allocatable, dimension(:) :: &
rho_mob_0, & !< initial dislocation density
rho_dip_0, & !< initial dipole density
b_sl, & !< magnitude of burgers vector [m]
@ -46,30 +35,30 @@ submodule(constitutive) plastic_disloUCLA
kink_height, & !< height of the kink pair
w, & !< width of the kink pair
omega !< attempt frequency for kink pair nucleation
real(pReal), dimension(:,:), allocatable :: &
real(pReal), allocatable, dimension(:,:) :: &
h_sl_sl, & !< slip resistance from slip activity
forestProjectionEdge
real(pReal), dimension(:,:,:), allocatable :: &
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid, &
nonSchmid_pos, &
nonSchmid_neg
integer :: &
sum_N_sl !< total number of active slip system
integer, dimension(:), allocatable :: &
integer, allocatable, dimension(:) :: &
N_sl !< number of active slip systems for each family
integer(kind(undefined_ID)), dimension(:),allocatable :: &
outputID !< ID of each post result output
character(len=pStringLen), allocatable, dimension(:) :: &
output
logical :: &
dipoleFormation !< flag indicating consideration of dipole formation
end type !< container type for internal constitutive parameters
type :: tDisloUCLAState
real(pReal), dimension(:,:), pointer :: &
rho_mob, &
rho_dip, &
gamma_sl
end type tDisloUCLAState
type :: tDisloUCLAdependentState
real(pReal), dimension(:,:), allocatable :: &
Lambda_sl, &
@ -88,41 +77,36 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module subroutine plastic_disloUCLA_init
integer :: &
Ninstance, &
p, i, &
NipcMyPhase, &
sizeState, sizeDotState, &
startIndex, endIndex
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_DISLOUCLA_label//' init -+>>>'
write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
write(6,'(a)') ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_DISLOUCLA_label//' init -+>>>'; flush(6)
write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
write(6,'(a)') ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
Ninstance = count(phase_plasticity == PLASTICITY_DISLOUCLA_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
allocate(dependentState(Ninstance))
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
@ -134,9 +118,9 @@ module subroutine plastic_disloUCLA_init
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%mu = lattice_mu(p)
prm%aTol_rho = config%getFloat('atol_rho')
! sanity checks
if (prm%aTol_rho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho'
@ -147,7 +131,7 @@ module subroutine plastic_disloUCLA_init
slipActive: if (prm%sum_N_sl > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
@ -158,20 +142,20 @@ module subroutine plastic_disloUCLA_init
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
endif
prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%forestProjectionEdge = lattice_forestProjection_edge(prm%N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%forestProjectionEdge = transpose(prm%forestProjectionEdge)
prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl))
prm%rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%N_sl))
prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl))
prm%b_sl = config%getFloats('slipburgers', requiredSize=size(prm%N_sl))
prm%delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl))
prm%i_sl = config%getFloats('clambdaslip', requiredSize=size(prm%N_sl))
prm%tau_0 = config%getFloats('tau_peierls', requiredSize=size(prm%N_sl))
prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl), &
@ -182,14 +166,14 @@ module subroutine plastic_disloUCLA_init
prm%w = config%getFloats('kink_width', requiredSize=size(prm%N_sl))
prm%omega = config%getFloats('omega', requiredSize=size(prm%N_sl))
prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl))
prm%D = config%getFloat('grainsize')
prm%D_0 = config%getFloat('d0')
prm%Q_cl = config%getFloat('qsd')
prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal
prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl
prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key
! expand: family => system
prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl)
prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl)
@ -206,8 +190,8 @@ module subroutine plastic_disloUCLA_init
prm%i_sl = math_expand(prm%i_sl, prm%N_sl)
prm%atomicVolume = math_expand(prm%atomicVolume, prm%N_sl)
prm%D_a = math_expand(prm%D_a, prm%N_sl)
! sanity checks
if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0'
if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
@ -219,7 +203,7 @@ module subroutine plastic_disloUCLA_init
if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0'
if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl'
if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipb_sl'
else slipActive
allocate(prm%rho_mob_0(0))
allocate(prm%rho_dip_0(0))
@ -232,39 +216,14 @@ module subroutine plastic_disloUCLA_init
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(trim(outputs(i)))
case ('edge_density')
outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl>0)
case ('dipole_density')
outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl>0)
case ('shear_rate','shearrate','shear_rate_slip','shearrate_slip')
outputID = merge(dot_gamma_sl_ID,undefined_ID,prm%sum_N_sl>0)
case ('accumulated_shear','accumulatedshear','accumulated_shear_slip')
outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl>0)
case ('mfp','mfp_slip')
outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl>0)
case ('threshold_stress','threshold_stress_slip')
outputID = merge(tau_pass_ID,undefined_ID,prm%sum_N_sl>0)
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID, outputID]
endif
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
@ -275,14 +234,14 @@ module subroutine plastic_disloUCLA_init
stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase)
dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase)
dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:)
@ -290,21 +249,21 @@ module subroutine plastic_disloUCLA_init
plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal ! Don't use for convergence check
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
allocate(dst%Lambda_sl(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal)
allocate(dst%threshold_stress(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_disloUCLA_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @brief Calculate plastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
Mp,T,instance,of)
@ -312,7 +271,7 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
@ -320,18 +279,18 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
integer, intent(in) :: &
instance, &
of
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%sum_N_sl) :: &
dot_gamma_pos,dot_gamma_neg, &
ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics(Mp,T,instance,of,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
do i = 1, prm%sum_N_sl
Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%Schmid(1:3,1:3,i)
@ -340,17 +299,17 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
+ ddot_gamma_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ ddot_gamma_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo
end associate
end subroutine plastic_disloUCLA_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
@ -358,7 +317,7 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
VacancyDiffusion
real(pReal), dimension(param(instance)%sum_N_sl) :: &
@ -369,16 +328,16 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of)
dot_rho_dip_formation, &
dot_rho_dip_climb, &
dip_distance
associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance))
call kinetics(Mp,T,instance,of,&
gdot_pos,gdot_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg)
dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs
VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg
dot_rho_dip_formation = 0.0_pReal
dot_rho_dip_climb = 0.0_pReal
@ -393,73 +352,73 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of)
* (1.0_pReal/(dip_distance+prm%D_a))
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency?
end where
dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%Lambda_sl(:,of)) & ! multiplication
- dot_rho_dip_formation &
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) ! Spontaneous annihilation of 2 single edge dislocations
dot%rho_dip(:,of) = dot_rho_dip_formation &
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent
- dot_rho_dip_climb
end associate
end subroutine plastic_disloUCLA_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!> @brief Calculate derived quantities from state.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_disloUCLA_dependentState(instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%sum_N_sl) :: &
dislocationSpacing
associate(prm => param(instance), stt => state(instance),dst => dependentState(instance))
dislocationSpacing = sqrt(matmul(prm%forestProjectionEdge, &
stt%rho_mob(:,of)+stt%rho_dip(:,of)))
dst%threshold_stress(:,of) = prm%mu*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of)))
dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
end associate
end subroutine plastic_disloUCLA_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_disloUCLA_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (rho_mob_ID)
call results_writeDataset(group,stt%rho_mob,'rho_mob',&
'mobile dislocation density','1/m²')
case (rho_dip_ID)
call results_writeDataset(group,stt%rho_dip,'rho_dip',&
'dislocation dipole density''1/m²')
case (dot_gamma_sl_ID)
call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& ! this is not dot!!
'plastic shear','1')
case (Lambda_sl_ID)
call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
'mean free path for slip','m')
case (tau_pass_ID)
call results_writeDataset(group,dst%threshold_stress,'tau_pass',&
'threshold stress for slip','Pa')
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('edge_density') ! ToDo: should be rho_mob
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,'rho_mob',&
'mobile dislocation density','1/m²')
case('dipole_density') ! ToDo: should be rho_dip
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',&
'dislocation dipole density''1/m²')
case('shear_rate_slip') ! should be gamma
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& ! this is not dot!!
'plastic shear','1')
case('mfp_slip') !ToDo: should be Lambda
if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
'mean free path for slip','m')
case('threshold_stress_slip') !ToDo: should be tau_pass
if(prm%sum_N_sl>0) call results_writeDataset(group,dst%threshold_stress,'tau_pass',&
'threshold stress for slip','Pa')
end select
enddo outputsLoop
end associate
@ -468,15 +427,15 @@ end subroutine plastic_disloUCLA_results
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the
! resolved stresss
!> @brief Calculate shear rates on slip systems, their derivatives with respect to resolved
! stress, and the resolved stress.
!> @details Derivatives and resolved stress are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,T,instance,of, &
dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
@ -484,7 +443,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
integer, intent(in) :: &
instance, &
of
real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
dot_gamma_pos, &
dot_gamma_neg
@ -501,82 +460,82 @@ pure subroutine kinetics(Mp,T,instance,of, &
t_n, t_k, dtk,dtn, &
needsGoodName ! ToDo: @Karo: any idea?
integer :: j
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
do j = 1, prm%sum_N_sl
tau_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j))
tau_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j))
enddo
if (present(tau_pos_out)) tau_pos_out = tau_pos
if (present(tau_neg_out)) tau_neg_out = tau_neg
associate(BoltzmannRatio => prm%delta_F/(kB*T), &
dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, &
effectiveLength => dst%Lambda_sl(:,of) - prm%w)
significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0
StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
vel = prm%kink_height/(t_n + t_k)
dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal
else where significantPositiveTau
dot_gamma_pos = 0.0_pReal
end where significantPositiveTau
if (present(ddot_gamma_dtau_pos)) then
significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0
dtk = -1.0_pReal * t_k / tau_pos
dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal
ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal
else where significantPositiveTau2
ddot_gamma_dtau_pos = 0.0_pReal
end where significantPositiveTau2
endif
significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_0
StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
vel = prm%kink_height/(t_n + t_k)
dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal
else where significantNegativeTau
dot_gamma_neg = 0.0_pReal
end where significantNegativeTau
if (present(ddot_gamma_dtau_neg)) then
significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0
dtk = -1.0_pReal * t_k / tau_neg
dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal
ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal
else where significantNegativeTau2
ddot_gamma_dtau_neg = 0.0_pReal
end where significantNegativeTau2
end if
end associate
end associate

519
src/constitutive_plastic_dislotwin.f90
View File

File diff suppressed because it is too large Load Diff

153
src/constitutive_plastic_isotropic.f90
View File

@ -8,14 +8,7 @@
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_isotropic
enum, bind(c)
enumerator :: &
undefined_ID, &
xi_ID, &
dot_gamma_ID
end enum
type :: tParameters
real(pReal) :: &
M, & !< Taylor factor
@ -34,12 +27,12 @@ submodule(constitutive) plastic_isotropic
aTol_gamma
integer :: &
of_debug = 0
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
logical :: &
dilatation
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type :: tIsotropicState
real(pReal), pointer, dimension(:) :: &
xi, &
@ -56,50 +49,45 @@ submodule(constitutive) plastic_isotropic
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_init
integer :: &
Ninstance, &
p, i, &
p, &
NipcMyPhase, &
sizeState, sizeDotState
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'
write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'; flush(6)
write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
Ninstance = count(phase_plasticity == PLASTICITY_ISOTROPIC_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), &
config => config_phase(p))
#ifdef DEBUG
if (p==material_phaseAt(debug_g,debug_e)) &
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
#endif
prm%xi_0 = config%getFloat('tau0')
prm%xi_inf = config%getFloat('tausat')
prm%dot_gamma_0 = config%getFloat('gdot0')
@ -114,9 +102,9 @@ module subroutine plastic_isotropic_init
prm%a = config%getFloat('a')
prm%aTol_xi = config%getFloat('atol_flowstress',defaultVal=1.0_pReal)
prm%aTol_gamma = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
prm%dilatation = config%keyExists('/dilatation/')
!--------------------------------------------------------------------------------------------------
! sanity checks
extmsg = ''
@ -128,79 +116,62 @@ module subroutine plastic_isotropic_init
if (prm%M <= 0.0_pReal) extmsg = trim(extmsg)//' m'
if (prm%aTol_xi <= 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
if (prm%aTol_gamma <= 0.0_pReal) extmsg = trim(extmsg)//' atol_shear'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('flowstress')
outputID = xi_ID
case ('strainrate')
outputID = dot_gamma_ID
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID, outputID]
endif
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['xi ','accumulated_shear'])
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
stt%xi => plasticState(p)%state (1,:)
stt%xi = prm%xi_0
dot%xi => plasticState(p)%dotState(1,:)
plasticState(p)%aTolState(1) = prm%aTol_xi
stt%gamma => plasticState(p)%state (2,:)
dot%gamma => plasticState(p)%dotState(2,:)
plasticState(p)%aTolState(2) = prm%aTol_gamma
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_isotropic_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @brief Calculate plastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(3,3) :: &
Mp_dev !< deviatoric part of the Mandel stress
real(pReal) :: &
@ -209,16 +180,16 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
squarenorm_Mp_dev !< square of the norm of the deviatoric part of the Mandel stress
integer :: &
k, l, m, n
associate(prm => param(instance), stt => state(instance))
Mp_dev = math_deviatoric33(Mp)
squarenorm_Mp_dev = math_mul33xx33(Mp_dev,Mp_dev)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
if (norm_Mp_dev > 0.0_pReal) then
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(of))) **prm%n
Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
@ -240,42 +211,42 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
end if
end associate
end subroutine plastic_isotropic_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @brief Calculate inelastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dMi !< derivative of Li with respect to Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mi !< Mandel stress
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
tr !< trace of spherical part of Mandel stress (= 3 x pressure)
integer :: &
k, l, m, n
associate(prm => param(instance), stt => state(instance))
tr=math_trace33(math_spherical33(Mi))
if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero
Li = math_I3 &
* prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) &
* tr * abs(tr)**(prm%n-1.0_pReal)
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
.and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then
@ -292,38 +263,38 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
Li = 0.0_pReal
dLi_dMi = 0.0_pReal
endif
end associate
end subroutine plastic_isotropic_LiAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
dot_gamma, & !< strainrate
xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
if (prm%dilatation) then
norm_Mp = sqrt(math_mul33xx33(Mp,Mp))
else
norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp)))
endif
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(of))) **prm%n
if (dot_gamma > 1e-12_pReal) then
if (dEq0(prm%c_1)) then
xi_inf_star = prm%xi_inf
@ -339,28 +310,28 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of)
else
dot%xi(of) = 0.0_pReal
endif
dot%gamma(of) = dot_gamma ! ToDo: not really used
end associate
end subroutine plastic_isotropic_dotState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_results(instance,group)
integer, intent(in) :: instance
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (xi_ID)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('flowstress') ! ToDo: should be 'xi'
call results_writeDataset(group,stt%xi,'xi','resistance against plastic flow','Pa')
end select
enddo outputsLoop

256
src/constitutive_plastic_kinehardening.f90
View File

@ -7,26 +7,13 @@
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_kinehardening
enum, bind(c)
enumerator :: &
undefined_ID, &
crss_ID, & !< critical resolved stress
crss_back_ID, & !< critical resolved back stress
sense_ID, & !< sense of acting shear stress (-1 or +1)
chi0_ID, & !< backstress at last switch of stress sense (positive?)
gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?)
accshear_ID, &
shearrate_ID, &
resolvedstress_ID
end enum
type :: tParameters
real(pReal) :: &
gdot0, & !< reference shear strain rate for slip
n, & !< stress exponent for slip
aTolResistance, &
aTolShear
real(pReal), allocatable, dimension(:) :: &
real(pReal), allocatable, dimension(:) :: &
crss0, & !< initial critical shear stress for slip
theta0, & !< initial hardening rate of forward stress for each slip
theta1, & !< asymptotic hardening rate of forward stress for each slip
@ -35,21 +22,21 @@ submodule(constitutive) plastic_kinehardening
tau1, &
tau1_b, &
nonSchmidCoeff
real(pReal), allocatable, dimension(:,:) :: &
real(pReal), allocatable, dimension(:,:) :: &
interaction_slipslip !< slip resistance from slip activity
real(pReal), allocatable, dimension(:,:,:) :: &
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid, &
nonSchmid_pos, &
nonSchmid_neg
integer :: &
totalNslip, & !< total number of active slip system
of_debug = 0
integer, allocatable, dimension(:) :: &
integer, allocatable, dimension(:) :: &
Nslip !< number of active slip systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress
@ -59,7 +46,7 @@ submodule(constitutive) plastic_kinehardening
gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear
end type tKinehardeningState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
@ -72,37 +59,32 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_init
integer :: &
Ninstance, &
p, i, o, &
p, o, &
NipcMyPhase, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'; flush(6)
Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
allocate(deltaState(Ninstance))
do p = 1, size(phase_plasticityInstance)
if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
@ -110,22 +92,22 @@ module subroutine plastic_kinehardening_init
dlt => deltaState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)),&
config => config_phase(p))
#ifdef DEBUG
if (p==material_phaseAt(debug_g,debug_e)) then
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
endif
#endif
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
! sanity checks
if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
@ -133,7 +115,7 @@ module subroutine plastic_kinehardening_init
slipActive: if (prm%totalNslip > 0) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
@ -146,7 +128,7 @@ module subroutine plastic_kinehardening_init
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip))
@ -154,10 +136,10 @@ module subroutine plastic_kinehardening_init
prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip))
prm%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip')
prm%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip')
! expand: family => system
prm%crss0 = math_expand(prm%crss0, prm%Nslip)
prm%tau1 = math_expand(prm%tau1, prm%Nslip)
@ -166,9 +148,9 @@ module subroutine plastic_kinehardening_init
prm%theta1 = math_expand(prm%theta1, prm%Nslip)
prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip)
prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip)
!--------------------------------------------------------------------------------------------------
! sanity checks
if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
@ -176,58 +158,29 @@ module subroutine plastic_kinehardening_init
if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0'
if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1'
if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b'
!ToDo: Any sensible checks for theta?
endif slipActive
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance')
outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0)
case ('accumulatedshear')
outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0)
case ('shearrate')
outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0)
case ('resolvedstress')
outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0)
case ('backstress')
outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0)
case ('sense')
outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0)
case ('chi0')
outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0)
case ('gamma0')
outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0)
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID , outputID]
endif
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip
sizeState = sizeDotState + sizeDeltaState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1
@ -236,13 +189,13 @@ module subroutine plastic_kinehardening_init
stt%crss = spread(prm%crss0, 2, NipcMyPhase)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%crss_back => plasticState(p)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%accshear => plasticState(p)%state (startIndex:endIndex,:)
@ -250,34 +203,34 @@ module subroutine plastic_kinehardening_init
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
o = plasticState(p)%offsetDeltaState
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @brief Calculate plastic velocity gradient and its tangent.
!--------------------------------------------------------------------------------------------------
pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
@ -285,26 +238,26 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1, prm%totalNslip
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -312,14 +265,14 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
+ dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo
end associate
end subroutine plastic_kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_dotState(Mp,instance,of)
@ -328,40 +281,40 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
sumGamma
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,of))
dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
* ( prm%theta1 &
+ (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) &
* exp(-sumGamma*prm%theta0/prm%tau1) &
)
dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
( prm%theta1_b + &
(prm%theta0_b - prm%theta1_b &
+ prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) &
)
end associate
end subroutine plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calculates (instantaneous) incremental change of microstructure
!> @brief Calculate (instantaneous) incremental change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
@ -370,18 +323,18 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_pos,gdot_neg, &
sense
associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
sense = merge(state(instance)%sense(:,of), & ! keep existing...
sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined
dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
.and. (of == prm%of_debug &
@ -390,7 +343,7 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
write(6,*) sense,state(instance)%sense(:,of)
endif
#endif
!--------------------------------------------------------------------------------------------------
! switch in sense of shear?
where(dNeq(sense,stt%sense(:,of),0.1_pReal))
@ -402,45 +355,43 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
dlt%chi0 (:,of) = 0.0_pReal
dlt%gamma0(:,of) = 0.0_pReal
end where
end associate
end subroutine plastic_kinehardening_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (crss_ID)
call results_writeDataset(group,stt%crss,'xi_sl', &
'resistance against plastic slip','Pa')
case(crss_back_ID)
call results_writeDataset(group,stt%crss_back,'tau_back', &
'back stress against plastic slip','Pa')
case (sense_ID)
call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1')
case (chi0_ID)
call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa')
case (gamma0_ID)
call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1')
case (accshear_ID)
call results_writeDataset(group,stt%accshear,'gamma_sl', &
'plastic shear','1')
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('resistance')
if(prm%totalNslip>0) call results_writeDataset(group,stt%crss,'xi_sl', &
'resistance against plastic slip','Pa')
case('backstress') ! ToDo: should be 'tau_back'
if(prm%totalNslip>0) call results_writeDataset(group,stt%crss_back,'tau_back', &
'back stress against plastic slip','Pa')
case ('sense')
if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear', &
'tbd','1')
case ('chi0')
if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0', &
'tbd','Pa')
case ('gamma0')
if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0', &
'tbd','1')
case ('accumulatedshear')
if(prm%totalNslip>0) call results_writeDataset(group,stt%accshear,'gamma_sl', &
'plastic shear','1')
end select
enddo outputsLoop
end associate
@ -449,10 +400,11 @@ end subroutine plastic_kinehardening_results
!--------------------------------------------------------------------------------------------------
!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress
!> @details: Shear rates are calculated only optionally.
!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved
! stress.
!> @details: Derivatives are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,instance,of, &
gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
@ -462,44 +414,44 @@ pure subroutine kinetics(Mp,instance,of, &
integer, intent(in) :: &
instance, &
of
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: &
gdot_pos, &
gdot_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: &
dgdot_dtau_pos, &
dgdot_dtau_neg
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_pos, &
tau_neg
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip
tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of)
tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), &
0.0_pReal, nonSchmidActive)
enddo
where(dNeq0(tau_pos))
gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos)
else where
gdot_pos = 0.0_pReal
end where
where(dNeq0(tau_neg))
gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg)
else where
gdot_neg = 0.0_pReal
end where
if (present(dgdot_dtau_pos)) then
where(dNeq0(gdot_pos))
dgdot_dtau_pos = gdot_pos*prm%n/tau_pos

12
src/constitutive_plastic_none.f90
View File

@ -18,19 +18,19 @@ module subroutine plastic_none_init
Ninstance, &
p, &
NipcMyPhase
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>'; flush(6)
Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
call material_allocatePlasticState(p,NipcMyPhase,0,0,0)
enddo
end subroutine plastic_none_init

203
src/constitutive_plastic_nonlocal.f90
View File

@ -48,32 +48,6 @@ submodule(constitutive) plastic_nonlocal
real(pReal), dimension(:,:,:,:,:,:), allocatable :: &
compatibility !< slip system compatibility between me and my neighbors
enum, bind(c)
enumerator :: &
undefined_ID, &
rho_sgl_mob_edg_pos_ID, &
rho_sgl_mob_edg_neg_ID, &
rho_sgl_mob_scr_pos_ID, &
rho_sgl_mob_scr_neg_ID, &
rho_sgl_imm_edg_pos_ID, &
rho_sgl_imm_edg_neg_ID, &
rho_sgl_imm_scr_pos_ID, &
rho_sgl_imm_scr_neg_ID, &
rho_dip_edg_ID, &
rho_dip_scr_ID, &
rho_forest_ID, &
resolvedstress_back_ID, &
tau_pass_ID, &
rho_dot_sgl_ID, &
rho_dot_sgl_mobile_ID, &
rho_dot_dip_ID, &
v_edg_pos_ID, &
v_edg_neg_ID, &
v_scr_pos_ID, &
v_scr_neg_ID, &
gamma_ID
end enum
type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
atomicVolume, & !< atomic volume
@ -135,14 +109,12 @@ submodule(constitutive) plastic_nonlocal
integer, dimension(:) ,allocatable :: &
Nslip,&
colinearSystem !< colinear system to the active slip system (only valid for fcc!)
character(len=pStringLen), allocatable, dimension(:) :: &
output
logical :: &
shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term
probabilisticMultiplication
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID !< ID of each post result output
end type tParameters
type :: tNonlocalMicrostructure
@ -198,22 +170,19 @@ module subroutine plastic_nonlocal_init
integer :: &
sizeState, sizeDotState,sizeDependentState, sizeDeltaState, &
maxNinstances, &
p, i, &
p, &
l, &
s1, s2, &
s, &
t, &
c
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = '', &
structure
character(len=pStringLen), dimension(:), allocatable :: outputs
integer :: NofMyPhase
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>'; flush(6)
write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014'
write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2014.03.012'
@ -407,60 +376,7 @@ module subroutine plastic_nonlocal_init
endif slipActive
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(trim(outputs(i)))
case ('rho_sgl_mob_edg_pos')
outputID = merge(rho_sgl_mob_edg_pos_ID,undefined_ID,prm%totalNslip>0)
case ('rho_sgl_mob_edg_neg')
outputID = merge(rho_sgl_mob_edg_neg_ID,undefined_ID,prm%totalNslip>0)
case ('rho_sgl_mob_scr_pos')
outputID = merge(rho_sgl_mob_scr_pos_ID,undefined_ID,prm%totalNslip>0)
case ('rho_sgl_mob_scr_neg')
outputID = merge(rho_sgl_mob_scr_neg_ID,undefined_ID,prm%totalNslip>0)
case ('rho_sgl_imm_edg_pos')
outputID = merge(rho_sgl_imm_edg_pos_ID,undefined_ID,prm%totalNslip>0)
case ('rho_sgl_imm_edg_neg')
outputID = merge(rho_sgl_imm_edg_neg_ID,undefined_ID,prm%totalNslip>0)
case ('rho_sgl_imm_scr_pos')
outputID = merge(rho_sgl_imm_scr_pos_ID,undefined_ID,prm%totalNslip>0)
case ('rho_sgl_imm_scr_neg')
outputID = merge(rho_sgl_imm_scr_neg_ID,undefined_ID,prm%totalNslip>0)
case ('rho_dip_edg')
outputID = merge(rho_dip_edg_ID,undefined_ID,prm%totalNslip>0)
case ('rho_dip_scr')
outputID = merge(rho_dip_scr_ID,undefined_ID,prm%totalNslip>0)
case ('rho_forest')
outputID = merge(rho_forest_ID,undefined_ID,prm%totalNslip>0)
case ('resolvedstress_back')
outputID = merge(resolvedstress_back_ID,undefined_ID,prm%totalNslip>0)
case ('tau_pass')
outputID = merge(tau_pass_ID,undefined_ID,prm%totalNslip>0)
case ('rho_dot_sgl')
outputID = merge(rho_dot_sgl_ID,undefined_ID,prm%totalNslip>0)
case ('rho_dot_sgl_mobile')
outputID = merge(rho_dot_sgl_mobile_ID,undefined_ID,prm%totalNslip>0)
case ('rho_dot_dip')
outputID = merge(rho_dot_dip_ID,undefined_ID,prm%totalNslip>0)
case ('v_edg_pos')
outputID = merge(v_edg_pos_ID,undefined_ID,prm%totalNslip>0)
case ('v_edg_neg')
outputID = merge(v_edg_neg_ID,undefined_ID,prm%totalNslip>0)
case ('v_scr_pos')
outputID = merge(v_scr_pos_ID,undefined_ID,prm%totalNslip>0)
case ('v_scr_neg')
outputID = merge(v_scr_neg_ID,undefined_ID,prm%totalNslip>0)
case ('gamma')
outputID = merge(gamma_ID,undefined_ID,prm%totalNslip>0)
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID , outputID]
endif
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
@ -561,7 +477,7 @@ module subroutine plastic_nonlocal_init
stt%v_scr_neg => plasticState(p)%state (15*prm%totalNslip + 1:16*prm%totalNslip ,1:NofMyPhase)
allocate(dst%tau_pass(prm%totalNslip,NofMyPhase),source=0.0_pReal)
allocate(dst%tau_Back(prm%totalNslip,NofMyPhase), source=0.0_pReal)
allocate(dst%tau_back(prm%totalNslip,NofMyPhase),source=0.0_pReal)
end associate
if (NofMyPhase > 0) call stateInit(p,NofMyPhase)
@ -1968,62 +1884,63 @@ module subroutine plastic_nonlocal_results(instance,group)
integer, intent(in) :: instance
character(len=*),intent(in) :: group
integer :: o
associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (rho_sgl_mob_edg_pos_ID)
call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', &
'positive mobile edge density','1/m²')
case (rho_sgl_imm_edg_pos_ID)
call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',&
'positive immobile edge density','1/m²')
case (rho_sgl_mob_edg_neg_ID)
call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',&
'negative mobile edge density','1/m²')
case (rho_sgl_imm_edg_neg_ID)
call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',&
'negative immobile edge density','1/m²')
case (rho_dip_edg_ID)
call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',&
'edge dipole density','1/m²')
case (rho_sgl_mob_scr_pos_ID)
call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',&
'positive mobile screw density','1/m²')
case (rho_sgl_imm_scr_pos_ID)
call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',&
'positive immobile screw density','1/m²')
case (rho_sgl_mob_scr_neg_ID)
call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',&
'negative mobile screw density','1/m²')
case (rho_sgl_imm_scr_neg_ID)
call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',&
'negative immobile screw density','1/m²')
case (rho_dip_scr_ID)
call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',&
'screw dipole density','1/m²')
case (rho_forest_ID)
call results_writeDataset(group,stt%rho_forest, 'rho_forest',&
'forest density','1/m²')
case (v_edg_pos_ID)
call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',&
'positive edge velocity','m/s')
case (v_edg_neg_ID)
call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',&
'negative edge velocity','m/s')
case (v_scr_pos_ID)
call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',&
'positive srew velocity','m/s')
case (v_scr_neg_ID)
call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',&
'negative screw velocity','m/s')
case(gamma_ID)
call results_writeDataset(group,stt%gamma,'gamma',&
'plastic shear','1')
case (tau_pass_ID)
call results_writeDataset(group,dst%tau_pass,'tau_pass',&
'passing stress for slip','Pa')
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('rho_sgl_mob_edg_pos')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', &
'positive mobile edge density','1/m²')
case('rho_sgl_imm_edg_pos')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',&
'positive immobile edge density','1/m²')
case('rho_sgl_mob_edg_neg')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',&
'negative mobile edge density','1/m²')
case('rho_sgl_imm_edg_neg')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',&
'negative immobile edge density','1/m²')
case('rho_dip_edg')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',&
'edge dipole density','1/m²')
case('rho_sgl_mob_scr_pos')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',&
'positive mobile screw density','1/m²')
case('rho_sgl_imm_scr_pos')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',&
'positive immobile screw density','1/m²')
case('rho_sgl_mob_scr_neg')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',&
'negative mobile screw density','1/m²')
case('rho_sgl_imm_scr_neg')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',&
'negative immobile screw density','1/m²')
case('rho_dip_scr')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',&
'screw dipole density','1/m²')
case('rho_forest')
if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_forest, 'rho_forest',&
'forest density','1/m²')
case('v_edg_pos')
if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',&
'positive edge velocity','m/s')
case('v_edg_neg')
if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',&
'negative edge velocity','m/s')
case('v_scr_pos')
if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',&
'positive srew velocity','m/s')
case('v_scr_neg')
if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',&
'negative screw velocity','m/s')
case('gamma')
if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma,'gamma',&
'plastic shear','1')
case('tau_pass')
if(prm%totalNslip>0) call results_writeDataset(group,dst%tau_pass,'tau_pass',&
'passing stress for slip','Pa')
end select
enddo outputsLoop
end associate

249
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -6,19 +6,6 @@
!--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_phenopowerlaw
enum, bind(c)
enumerator :: &
undefined_ID, &
resistance_slip_ID, &
accumulatedshear_slip_ID, &
shearrate_slip_ID, &
resolvedstress_slip_ID, &
resistance_twin_ID, &
accumulatedshear_twin_ID, &
shearrate_twin_ID, &
resolvedstress_twin_ID
end enum
type :: tParameters
real(pReal) :: &
gdot0_slip, & !< reference shear strain rate for slip
@ -37,19 +24,19 @@ submodule(constitutive) plastic_phenopowerlaw
aTolResistance, & !< absolute tolerance for integration of xi
aTolShear, & !< absolute tolerance for integration of gamma
aTolTwinfrac !< absolute tolerance for integration of f
real(pReal), allocatable, dimension(:) :: &
real(pReal), allocatable, dimension(:) :: &
xi_slip_0, & !< initial critical shear stress for slip
xi_twin_0, & !< initial critical shear stress for twin
xi_slip_sat, & !< maximum critical shear stress for slip
nonSchmidCoeff, &
H_int, & !< per family hardening activity (optional)
gamma_twin_char !< characteristic shear for twins
real(pReal), allocatable, dimension(:,:) :: &
real(pReal), allocatable, dimension(:,:) :: &
interaction_SlipSlip, & !< slip resistance from slip activity
interaction_SlipTwin, & !< slip resistance from twin activity
interaction_TwinSlip, & !< twin resistance from slip activity
interaction_TwinTwin !< twin resistance from twin activity
real(pReal), allocatable, dimension(:,:,:) :: &
real(pReal), allocatable, dimension(:,:,:) :: &
Schmid_slip, &
Schmid_twin, &
nonSchmid_pos, &
@ -57,13 +44,13 @@ submodule(constitutive) plastic_phenopowerlaw
integer :: &
totalNslip, & !< total number of active slip system
totalNtwin !< total number of active twin systems
integer, allocatable, dimension(:) :: &
integer, allocatable, dimension(:) :: &
Nslip, & !< number of active slip systems for each family
Ntwin !< number of active twin systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
xi_slip, &
@ -71,7 +58,7 @@ submodule(constitutive) plastic_phenopowerlaw
gamma_slip, &
gamma_twin
end type tPhenopowerlawState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
@ -83,7 +70,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Perform module initialization.
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_init
@ -91,51 +78,46 @@ module subroutine plastic_phenopowerlaw_init
integer :: &
Ninstance, &
p, i, &
NipcMyPhase, outputSize, &
NipcMyPhase, &
sizeState, sizeDotState, &
startIndex, endIndex
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'; flush(6)
Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(param(Ninstance))
allocate(state(Ninstance))
allocate(dotState(Ninstance))
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), &
config => config_phase(p))
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal)
prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal)
prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal)
prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal)
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal)
! sanity checks
if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac'
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
@ -143,7 +125,7 @@ module subroutine plastic_phenopowerlaw_init
slipActive: if (prm%totalNslip > 0) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
@ -157,22 +139,22 @@ module subroutine plastic_phenopowerlaw_init
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
config%getString('lattice_structure'))
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))])
prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip')
prm%a_slip = config%getFloat('a_slip')
prm%h0_SlipSlip = config%getFloat('h0_slipslip')
! expand: family => system
prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip)
prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip)
prm%H_int = math_expand(prm%H_int, prm%Nslip)
! sanity checks
if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip'
if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip'
@ -183,7 +165,7 @@ module subroutine plastic_phenopowerlaw_init
allocate(prm%interaction_SlipSlip(0,0))
allocate(prm%xi_slip_0(0))
endif slipActive
!--------------------------------------------------------------------------------------------------
! twin related parameters
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
@ -196,17 +178,17 @@ module subroutine plastic_phenopowerlaw_init
config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
prm%gdot0_twin = config%getFloat('gdot0_twin')
prm%n_twin = config%getFloat('n_twin')
prm%spr = config%getFloat('s_pr')
prm%h0_TwinTwin = config%getFloat('h0_twintwin')
! expand: family => system
prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin)
! sanity checks
if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin'
if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin'
@ -215,7 +197,7 @@ module subroutine plastic_phenopowerlaw_init
allocate(prm%xi_twin_0(0))
allocate(prm%gamma_twin_char(0))
endif twinActive
!--------------------------------------------------------------------------------------------------
! slip-twin related parameters
slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then
@ -231,63 +213,25 @@ module subroutine plastic_phenopowerlaw_init
allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0
prm%h0_TwinSlip = 0.0_pReal
endif slipAndTwinActive
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('resistance_slip')
outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('accumulatedshear_slip')
outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('shearrate_slip')
outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resolvedstress_slip')
outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0)
outputSize = prm%totalNslip
case ('resistance_twin')
outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('accumulatedshear_twin')
outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('shearrate_twin')
outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
case ('resolvedstress_twin')
outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0)
outputSize = prm%totalNtwin
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID, outputID]
endif
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip &
+ size(['tau_twin ','gamma_twin']) * prm%totalNtwin
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1
@ -296,14 +240,14 @@ module subroutine plastic_phenopowerlaw_init
stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNtwin
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1
endIndex = endIndex + prm%totalNslip
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
@ -317,18 +261,18 @@ module subroutine plastic_phenopowerlaw_init
stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
enddo
end subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @brief Calculate plastic velocity gradient and its tangent.
!> @details asummes that deformation by dislocation glide affects twinned and untwinned volume
! equally (Taylor assumption). Twinning happens only in untwinned volume
!--------------------------------------------------------------------------------------------------
@ -338,13 +282,13 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
integer :: &
i,k,l,m,n
real(pReal), dimension(param(instance)%totalNslip) :: &
@ -352,12 +296,12 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
dgdot_dtauslip_pos,dgdot_dtauslip_neg
real(pReal), dimension(param(instance)%totalNtwin) :: &
gdot_twin,dgdot_dtautwin
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1, prm%totalNslip
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
@ -366,7 +310,7 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
+ dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems
call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1, prm%totalNtwin
Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i)
@ -374,14 +318,14 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i)
enddo twinSystems
end associate
end subroutine plastic_phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
@ -390,7 +334,7 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal) :: &
c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,&
@ -398,9 +342,9 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
real(pReal), dimension(param(instance)%totalNslip) :: &
left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
sumGamma = sum(stt%gamma_slip(:,of))
sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)
@ -409,26 +353,26 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2)
c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3
c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4
!--------------------------------------------------------------------------------------------------
! calculate left and right vectors
left_SlipSlip = 1.0_pReal + prm%H_int
xi_slip_sat_offset = prm%spr*sqrt(sumF)
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset))
!--------------------------------------------------------------------------------------------------
! shear rates
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg)
dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg)
call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of))
!--------------------------------------------------------------------------------------------------
! hardening
dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * &
matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) &
+ matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of))
dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) &
+ c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of))
end associate
@ -437,33 +381,33 @@ end subroutine plastic_phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!> @brief Write results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_phenopowerlaw_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(instance), stt => state(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('resistance_slip')
if(prm%totalNslip>0) call results_writeDataset(group,stt%xi_slip, 'xi_sl', &
'resistance against plastic slip','Pa')
case('accumulatedshear_slip')
if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
'plastic shear','1')
case('resistance_twin')
if(prm%totalNtwin>0) call results_writeDataset(group,stt%xi_twin, 'xi_tw', &
'resistance against twinning','Pa')
case('accumulatedshear_twin')
if(prm%totalNtwin>0) call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
'twinning shear','1')
case (resistance_slip_ID)
call results_writeDataset(group,stt%xi_slip, 'xi_sl', &
'resistance against plastic slip','Pa')
case (accumulatedshear_slip_ID)
call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
'plastic shear','1')
case (resistance_twin_ID)
call results_writeDataset(group,stt%xi_twin, 'xi_tw', &
'resistance against twinning','Pa')
case (accumulatedshear_twin_ID)
call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
'twinning shear','1')
end select
enddo outputsLoop
end associate
@ -472,10 +416,11 @@ end subroutine plastic_phenopowerlaw_results
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress
!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved
! stress.
!> @details Derivatives are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,instance,of, &
gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
@ -485,44 +430,44 @@ pure subroutine kinetics_slip(Mp,instance,of, &
integer, intent(in) :: &
instance, &
of
real(pReal), intent(out), dimension(param(instance)%totalNslip) :: &
gdot_slip_pos, &
gdot_slip_neg
real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: &
dgdot_dtau_slip_pos, &
dgdot_dtau_slip_neg
real(pReal), dimension(param(instance)%totalNslip) :: &
tau_slip_pos, &
tau_slip_neg
integer :: i
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0
do i = 1, prm%totalNslip
tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
0.0_pReal, nonSchmidActive)
enddo
where(dNeq0(tau_slip_pos))
gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos)
else where
gdot_slip_pos = 0.0_pReal
end where
where(dNeq0(tau_slip_neg))
gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg)
else where
gdot_slip_neg = 0.0_pReal
end where
if (present(dgdot_dtau_slip_pos)) then
where(dNeq0(gdot_slip_pos))
dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos
@ -543,9 +488,9 @@ end subroutine kinetics_slip
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress.
! twinning is assumed to take place only in untwinned volume.
!> @details Derivates are calculated only optionally.
!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved
! stress. Twinning is assumed to take place only in untwinned volume.
!> @details Derivatives are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
@ -557,29 +502,29 @@ pure subroutine kinetics_twin(Mp,instance,of,&
integer, intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: &
gdot_twin
real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: &
dgdot_dtau_twin
real(pReal), dimension(param(instance)%totalNtwin) :: &
tau_twin
integer :: i
associate(prm => param(instance), stt => state(instance))
do i = 1, prm%totalNtwin
tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
enddo
where(tau_twin > 0.0_pReal)
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction
* prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
else where
gdot_twin = 0.0_pReal
end where
if (present(dgdot_dtau_twin)) then
where(dNeq0(gdot_twin))
dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin
@ -587,7 +532,7 @@ pure subroutine kinetics_twin(Mp,instance,of,&
dgdot_dtau_twin = 0.0_pReal
end where
endif
end associate
end subroutine kinetics_twin

70
src/crystallite.f90
View File

@ -57,9 +57,7 @@ module crystallite
crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc
real(pReal), dimension(:,:,:,:,:), allocatable :: &
crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step)
crystallite_subFp0,& !< plastic def grad at start of crystallite inc
crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step)
crystallite_subFi0,& !< intermediate def grad at start of crystallite inc
crystallite_subF, & !< def grad to be reached at end of crystallite inc
crystallite_subF0, & !< def grad at start of crystallite inc
@ -145,12 +143,10 @@ subroutine crystallite_init
allocate(crystallite_partionedFp0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subFp0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Fp(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_invFp(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Fi0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_partionedFi0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subFi0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Fi(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_invFi(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Fe(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Lp0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_partionedLp0(3,3,cMax,iMax,eMax), source=0.0_pReal)
@ -408,9 +404,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
else
crystallite_subStep(c,i,e) = num%subStepSizeCryst * crystallite_subStep(c,i,e)
crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e)
crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp (1:3,1:3,c,i,e))
crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e)
crystallite_invFi(1:3,1:3,c,i,e) = math_inv33(crystallite_Fi (1:3,1:3,c,i,e))
crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e)
if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback
crystallite_Lp (1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e)
@ -431,11 +425,11 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
! prepare for integration
if (crystallite_todo(c,i,e)) then
crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) &
+ crystallite_subStep(c,i,e) * (crystallite_partionedF (1:3,1:3,c,i,e) &
- crystallite_partionedF0(1:3,1:3,c,i,e))
crystallite_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e)), &
crystallite_invFi(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_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
math_inv33(crystallite_Fp(1:3,1:3,c,i,e))), &
math_inv33(crystallite_Fi(1:3,1:3,c,i,e)))
crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e)
crystallite_converged(c,i,e) = .false.
endif
@ -477,7 +471,9 @@ subroutine crystallite_stressTangent
o, &
p
real(pReal), dimension(3,3) :: temp_33_1, devNull,invSubFi0, temp_33_2, temp_33_3, temp_33_4
real(pReal), dimension(3,3) :: devNull, &
invSubFp0,invSubFi0,invFp,invFi, &
temp_33_1, temp_33_2, temp_33_3, temp_33_4
real(pReal), dimension(3,3,3,3) :: dSdFe, &
dSdF, &
dSdFi, &
@ -493,7 +489,8 @@ subroutine crystallite_stressTangent
real(pReal), dimension(9,9):: temp_99
logical :: error
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,invSubFi0,o,p, &
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,o,p, &
!$OMP invSubFp0,invSubFi0,invFp,invFi, &
!$OMP rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error)
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
@ -507,16 +504,20 @@ subroutine crystallite_stressTangent
crystallite_Fi(1:3,1:3,c,i,e), &
c,i,e)
invFp = math_inv33(crystallite_Fp(1:3,1:3,c,i,e))
invFi = math_inv33(crystallite_Fi(1:3,1:3,c,i,e))
invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))
invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))
if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal
else
invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))
lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
do o=1,3; do p=1,3
lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
+ crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
+ crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e)
+ invFi*invFi(p,o)
rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
- crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
enddo; enddo
@ -538,18 +539,13 @@ subroutine crystallite_stressTangent
!--------------------------------------------------------------------------------------------------
! calculate dSdF
temp_33_1 = transpose(matmul(crystallite_invFp(1:3,1:3,c,i,e), &
crystallite_invFi(1:3,1:3,c,i,e)))
temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e), &
math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)))
temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp (1:3,1:3,c,i,e)), &
math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)))
temp_33_1 = transpose(matmul(invFp,invFi))
temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e),invSubFp0)
temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),invFp), invSubFi0)
do o=1,3; do p=1,3
rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), &
crystallite_invFi(1:3,1:3,c,i,e)) &
temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
+ matmul(temp_33_3,dLidS(1:3,1:3,p,o))
enddo; enddo
lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) &
@ -569,15 +565,14 @@ subroutine crystallite_stressTangent
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
do o=1,3; do p=1,3
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e) &
* matmul(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), &
matmul(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e)))
* matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi))
enddo; enddo
!--------------------------------------------------------------------------------------------------
! assemble dPdF
temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(crystallite_invFp(1:3,1:3,c,i,e)))
temp_33_2 = matmul(crystallite_invFp(1:3,1:3,c,i,e),temp_33_1)
temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),crystallite_invFp(1:3,1:3,c,i,e))
temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(invFp))
temp_33_2 = matmul(invFp,temp_33_1)
temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),invFp)
temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,c,i,e))
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal
@ -589,7 +584,7 @@ subroutine crystallite_stressTangent
+ matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), &
transpose(crystallite_invFp(1:3,1:3,c,i,e))) &
transpose(invFp)) &
+ matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo
@ -793,28 +788,28 @@ logical function integrateStress(ipc,ip,el,timeFraction)
real(pReal), dimension(3,3):: F, & ! deformation gradient at end of timestep
Fp_new, & ! plastic deformation gradient at end of timestep
Fe_new, & ! elastic deformation gradient at end of timestep
invFp_new, & ! inverse of Fp_new
Fi_new, & ! gradient of intermediate deformation stages
invFi_new, &
invFp_current, & ! inverse of Fp_current
invFi_current, & ! inverse of Fp_current
Lpguess, & ! current guess for plastic velocity gradient
Lpguess_old, & ! known last good guess for plastic velocity gradient
Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law
residuumLp, & ! current residuum of plastic velocity gradient
residuumLp_old, & ! last residuum of plastic velocity gradient
deltaLp, & ! direction of next guess
Fi_new, & ! gradient of intermediate deformation stages
invFi_new, &
invFi_current, & ! inverse of Fi_current
Liguess, & ! current guess for intermediate velocity gradient
Liguess_old, & ! known last good guess for intermediate velocity gradient
Li_constitutive, & ! intermediate velocity gradient resulting from constitutive law
residuumLi, & ! current residuum of intermediate velocity gradient
residuumLi_old, & ! last residuum of intermediate velocity gradient
deltaLi, & ! direction of next guess
Fe, & ! elastic deformation gradient
Fe_new, &
S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
A, &
B, &
Fe, & ! elastic deformation gradient
temp_33
real(pReal), dimension(9) :: temp_9 ! needed for matrix inversion by LAPACK
integer, dimension(9) :: devNull_9 ! needed for matrix inversion by LAPACK
@ -993,16 +988,13 @@ logical function integrateStress(ipc,ip,el,timeFraction)
Fe_new = matmul(matmul(F,invFp_new),invFi_new)
integrateStress = .true.
crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) ! ToDo: We propably do not need to store P!
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_Lp (1:3,1:3,ipc,ip,el) = Lpguess
crystallite_Li (1:3,1:3,ipc,ip,el) = Liguess
crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new
crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new
crystallite_Fe (1:3,1:3,ipc,ip,el) = Fe_new
crystallite_invFp(1:3,1:3,ipc,ip,el) = invFp_new
crystallite_invFi(1:3,1:3,ipc,ip,el) = invFi_new
end function integrateStress

80
src/homogenization_mech_RGC.f90
View File

@ -9,17 +9,6 @@
submodule(homogenization) homogenization_mech_RGC
use rotations
enum, bind(c)
enumerator :: &
undefined_ID, &
constitutivework_ID, &
penaltyenergy_ID, &
volumediscrepancy_ID, &
averagerelaxrate_ID,&
maximumrelaxrate_ID,&
magnitudemismatch_ID
end enum
type :: tParameters
integer, dimension(:), allocatable :: &
Nconstituents
@ -31,8 +20,8 @@ submodule(homogenization) homogenization_mech_RGC
angles
integer :: &
of_debug = 0
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type :: tRGCstate
@ -71,23 +60,18 @@ module subroutine mech_RGC_init
integer :: &
Ninstance, &
h, i, &
h, &
NofMyHomog, &
sizeState, nIntFaceTot
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>'; flush(6)
write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
write(6,'(a)') ' https://doi.org/10.1007/s12289-009-0619-1'
write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
write(6,'(a)') ' https://doi.org/10.1007/s12289-009-0619-1'
write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
write(6,'(a)') ' https://doi.org/10.1088/0965-0393/18/1/015006'
write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
write(6,'(a)') ' https://doi.org/10.1088/0965-0393/18/1/015006'
Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID)
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) &
@ -123,34 +107,8 @@ module subroutine mech_RGC_init
prm%dAlpha = config%getFloats('grainsize', requiredSize=3)
prm%angles = config%getFloats('clusterorientation',requiredSize=3)
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case('constitutivework')
outputID = constitutivework_ID
case('penaltyenergy')
outputID = penaltyenergy_ID
case('volumediscrepancy')
outputID = volumediscrepancy_ID
case('averagerelaxrate')
outputID = averagerelaxrate_ID
case('maximumrelaxrate')
outputID = maximumrelaxrate_ID
case('magnitudemismatch')
outputID = magnitudemismatch_ID
end select
if (outputID /= undefined_ID) then
prm%outputID = [prm%outputID , outputID]
endif
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
NofMyHomog = count(material_homogenizationAt == h)
nIntFaceTot = 3*( (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) &
+ prm%Nconstituents(1)*(prm%Nconstituents(2)-1)*prm%Nconstituents(3) &
@ -934,26 +892,24 @@ module subroutine mech_RGC_results(instance,group)
integer :: o
associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (constitutivework_ID)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('constitutivework')
call results_writeDataset(group,stt%work,'W',&
'work density','J/m³')
case (magnitudemismatch_ID)
case('magnitudemismatch')
call results_writeDataset(group,dst%mismatch,'N',&
'average mismatch tensor','1')
case (penaltyenergy_ID)
case('penaltyenergy')
call results_writeDataset(group,stt%penaltyEnergy,'R',&
'mismatch penalty density','J/m³')
case (volumediscrepancy_ID)
case('volumediscrepancy')
call results_writeDataset(group,dst%volumeDiscrepancy,'Delta_V',&
'volume discrepancy','m³')
case (maximumrelaxrate_ID)
case('maximumrelaxrate')
call results_writeDataset(group,dst%relaxationrate_max,'max_alpha_dot',&
'maximum relaxation rate','m/s')
case (averagerelaxrate_ID)
case('averagerelaxrate')
call results_writeDataset(group,dst%relaxationrate_avg,'avg_alpha_dot',&
'average relaxation rate','m/s')
end select

34
src/thermal_adiabatic.f90
View File

@ -16,14 +16,9 @@ module thermal_adiabatic
implicit none
private
enum, bind(c)
enumerator :: undefined_ID, &
temperature_ID
end enum
type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tparameters), dimension(:), allocatable :: &
@ -46,10 +41,9 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_init
integer :: maxNinstance,o,h,NofMyHomog
character(len=pStringLen), dimension(:), allocatable :: outputs
integer :: maxNinstance,h,NofMyHomog
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>'; flush(6)
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>'; flush(6)
maxNinstance = count(thermal_type == THERMAL_adiabatic_ID)
if (maxNinstance == 0) return
@ -60,15 +54,7 @@ subroutine thermal_adiabatic_init
if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle
associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h))
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do o=1, size(outputs)
select case(outputs(o))
case('temperature')
prm%outputID = [prm%outputID, temperature_ID]
end select
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
NofMyHomog=count(material_homogenizationAt==h)
thermalState(h)%sizeState = 1
@ -76,7 +62,6 @@ subroutine thermal_adiabatic_init
allocate(thermalState(h)%subState0(1,NofMyHomog), source=thermal_initialT(h))
allocate(thermalState(h)%state (1,NofMyHomog), source=thermal_initialT(h))
nullify(thermalMapping(h)%p)
thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature(h)%p)
temperature(h)%p => thermalState(h)%state(1,:)
@ -246,14 +231,13 @@ subroutine thermal_adiabatic_results(homog,group)
integer, intent(in) :: homog
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(damage_typeInstance(homog)))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (temperature_ID)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('temperature') ! ToDo: should be 'T'
call results_writeDataset(group,temperature(homog)%p,'T',&
'temperature','K')
end select

35
src/thermal_conduction.f90
View File

@ -15,15 +15,9 @@ module thermal_conduction
implicit none
private
enum, bind(c)
enumerator :: &
undefined_ID, &
temperature_ID
end enum
type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tparameters), dimension(:), allocatable :: &
@ -48,10 +42,9 @@ contains
subroutine thermal_conduction_init
integer :: maxNinstance,o,NofMyHomog,h
character(len=pStringLen), dimension(:), allocatable :: outputs
integer :: maxNinstance,NofMyHomog,h
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'; flush(6)
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'; flush(6)
maxNinstance = count(thermal_type == THERMAL_conduction_ID)
if (maxNinstance == 0) return
@ -62,15 +55,7 @@ subroutine thermal_conduction_init
if (thermal_type(h) /= THERMAL_conduction_ID) cycle
associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h))
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do o=1, size(outputs)
select case(outputs(o))
case('temperature')
prm%outputID = [prm%outputID, temperature_ID]
end select
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
NofMyHomog=count(material_homogenizationAt==h)
thermalState(h)%sizeState = 0
@ -78,7 +63,6 @@ subroutine thermal_conduction_init
allocate(thermalState(h)%subState0(0,NofMyHomog))
allocate(thermalState(h)%state (0,NofMyHomog))
nullify(thermalMapping(h)%p)
thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature (h)%p)
allocate (temperature (h)%p(NofMyHomog), source=thermal_initialT(h))
@ -259,14 +243,13 @@ subroutine thermal_conduction_results(homog,group)
integer, intent(in) :: homog
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(damage_typeInstance(homog)))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (temperature_ID)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('temperature') ! ToDo: should be 'T'
call results_writeDataset(group,temperature(homog)%p,'T',&
'temperature','K')
end select