DAMASK with grid solvers  Revision: v2.0.3-2204-gdb1f2151
The Düsseldorf Advanced Material Simulation Kit with Grid Solvers
constitutive_plastic_phenopowerlaw.f90
Go to the documentation of this file.
1 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_phenopowerlaw.f90"
2 # 1 "<built-in>"
3 # 1 "<command-line>"
4 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_phenopowerlaw.f90"
5 !--------------------------------------------------------------------------------------------------
10 !--------------------------------------------------------------------------------------------------
11 submodule(constitutive) plastic_phenopowerlaw
12 
13  type :: tparameters
14  real(pReal) :: &
15  gdot0_slip = 1.0_preal, &
16  gdot0_twin = 1.0_preal, &
17  n_slip = 1.0_preal, &
18  n_twin = 1.0_preal, &
19  spr = 1.0_preal, &
20  c_1 = 1.0_preal, &
21  c_2 = 1.0_preal, &
22  c_3 = 1.0_preal, &
23  c_4 = 1.0_preal, &
24  h0_slipslip = 1.0_preal, &
25  h0_twinslip = 1.0_preal, &
26  h0_twintwin = 1.0_preal, &
27  a_slip = 1.0_preal
28  real(pReal), allocatable, dimension(:) :: &
29  xi_slip_sat, & !< maximum critical shear stress for slip
30  H_int, & !< per family hardening activity (optional)
31  gamma_twin_char
32  real(pReal), allocatable, dimension(:,:) :: &
33  interaction_SlipSlip, & !< slip resistance from slip activity
34  interaction_SlipTwin, & !< slip resistance from twin activity
35  interaction_TwinSlip, & !< twin resistance from slip activity
36  interaction_TwinTwin
37  real(pReal), allocatable, dimension(:,:,:) :: &
38  P_sl, &
39  P_tw, &
40  nonSchmid_pos, &
41  nonSchmid_neg
42  integer :: &
43  sum_N_sl, & !< total number of active slip system
44  sum_N_tw
45  logical :: &
46  nonSchmidActive = .false.
47  character(len=pStringLen), allocatable, dimension(:) :: &
48  output
49  end type tparameters
50 
51  type :: tphenopowerlawstate
52  real(pReal), pointer, dimension(:,:) :: &
53  xi_slip, &
54  xi_twin, &
55  gamma_slip, &
56  gamma_twin
57  end type tphenopowerlawstate
58 
59 !--------------------------------------------------------------------------------------------------
60 ! containers for parameters and state
61  type(tParameters), allocatable, dimension(:) :: param
62  type(tPhenopowerlawState), allocatable, dimension(:) :: &
63  dotState, &
64  state
65 
66 contains
67 
68 
69 !--------------------------------------------------------------------------------------------------
72 !--------------------------------------------------------------------------------------------------
73 module subroutine plastic_phenopowerlaw_init
74 
75  integer :: &
76  Ninstance, &
77  p, i, &
78  NipcMyPhase, &
79  sizeState, sizeDotState, &
80  startIndex, endIndex
81  integer, dimension(:), allocatable :: &
82  N_sl, N_tw
83  real(pReal), dimension(:), allocatable :: &
84  xi_slip_0, & !< initial critical shear stress for slip
85  xi_twin_0, & !< initial critical shear stress for twin
86  a
87  character(len=pStringLen) :: &
88  extmsg = ''
89 
90  write(6,'(/,a)') ' <<<+- plastic_'//plasticity_phenopowerlaw_label//' init -+>>>'; flush(6)
91 
92  ninstance = count(phase_plasticity == plasticity_phenopowerlaw_id)
93  if (iand(debug_level(debug_constitutive),debug_levelbasic) /= 0) &
94  write(6,'(a16,1x,i5,/)') '# instances:',ninstance
95 
96  allocate(param(ninstance))
97  allocate(state(ninstance))
98  allocate(dotstate(ninstance))
99 
100  do p = 1, size(phase_plasticity)
101  if (phase_plasticity(p) /= plasticity_phenopowerlaw_id) cycle
102  associate(prm => param(phase_plasticityinstance(p)), &
103  dot => dotstate(phase_plasticityinstance(p)), &
104  stt => state(phase_plasticityinstance(p)), &
105  config => config_phase(p))
106 
107 !--------------------------------------------------------------------------------------------------
108 ! slip related parameters
109  n_sl = config%getInts('nslip',defaultval=emptyintarray)
110  prm%sum_N_sl = sum(abs(n_sl))
111  slipactive: if (prm%sum_N_sl > 0) then
112  prm%P_sl = lattice_schmidmatrix_slip(n_sl,config%getString('lattice_structure'),&
113  config%getFloat('c/a',defaultval=0.0_preal))
114 
115  if(trim(config%getString('lattice_structure')) == 'bcc') then
116  a = config%getFloats('nonschmid_coefficients',defaultval = emptyrealarray)
117  if(size(a) > 0) prm%nonSchmidActive = .true.
118  prm%nonSchmid_pos = lattice_nonschmidmatrix(n_sl,a,+1)
119  prm%nonSchmid_neg = lattice_nonschmidmatrix(n_sl,a,-1)
120  else
121  prm%nonSchmid_pos = prm%P_sl
122  prm%nonSchmid_neg = prm%P_sl
123  endif
124  prm%interaction_SlipSlip = lattice_interaction_slipbyslip(n_sl, &
125  config%getFloats('interaction_slipslip'), &
126  config%getString('lattice_structure'))
127 
128  xi_slip_0 = config%getFloats('tau0_slip', requiredsize=size(n_sl))
129  prm%xi_slip_sat = config%getFloats('tausat_slip', requiredsize=size(n_sl))
130  prm%H_int = config%getFloats('h_int', requiredsize=size(n_sl), &
131  defaultval=[(0.0_preal,i=1,size(n_sl))])
132 
133  prm%gdot0_slip = config%getFloat('gdot0_slip')
134  prm%n_slip = config%getFloat('n_slip')
135  prm%a_slip = config%getFloat('a_slip')
136  prm%h0_SlipSlip = config%getFloat('h0_slipslip')
137 
138  ! expand: family => system
139  xi_slip_0 = math_expand(xi_slip_0, n_sl)
140  prm%xi_slip_sat = math_expand(prm%xi_slip_sat,n_sl)
141  prm%H_int = math_expand(prm%H_int, n_sl)
142 
143  ! sanity checks
144  if ( prm%gdot0_slip <= 0.0_preal) extmsg = trim(extmsg)//' gdot0_slip'
145  if ( prm%a_slip <= 0.0_preal) extmsg = trim(extmsg)//' a_slip'
146  if ( prm%n_slip <= 0.0_preal) extmsg = trim(extmsg)//' n_slip'
147  if (any(xi_slip_0 <= 0.0_preal)) extmsg = trim(extmsg)//' xi_slip_0'
148  if (any(prm%xi_slip_sat <= 0.0_preal)) extmsg = trim(extmsg)//' xi_slip_sat'
149 
150  else slipactive
151  xi_slip_0 = emptyrealarray
152  allocate(prm%xi_slip_sat,prm%H_int,source=emptyrealarray)
153  allocate(prm%interaction_SlipSlip(0,0))
154  endif slipactive
155 
156 !--------------------------------------------------------------------------------------------------
157 ! twin related parameters
158  n_tw = config%getInts('ntwin', defaultval=emptyintarray)
159  prm%sum_N_tw = sum(abs(n_tw))
160  twinactive: if (prm%sum_N_tw > 0) then
161  prm%P_tw = lattice_schmidmatrix_twin(n_tw,config%getString('lattice_structure'),&
162  config%getFloat('c/a',defaultval=0.0_preal))
163  prm%interaction_TwinTwin = lattice_interaction_twinbytwin(n_tw,&
164  config%getFloats('interaction_twintwin'), &
165  config%getString('lattice_structure'))
166  prm%gamma_twin_char = lattice_characteristicshear_twin(n_tw,config%getString('lattice_structure'),&
167  config%getFloat('c/a'))
168 
169  xi_twin_0 = config%getFloats('tau0_twin',requiredsize=size(n_tw))
170 
171  prm%c_1 = config%getFloat('twin_c',defaultval=0.0_preal)
172  prm%c_2 = config%getFloat('twin_b',defaultval=1.0_preal)
173  prm%c_3 = config%getFloat('twin_e',defaultval=0.0_preal)
174  prm%c_4 = config%getFloat('twin_d',defaultval=0.0_preal)
175  prm%gdot0_twin = config%getFloat('gdot0_twin')
176  prm%n_twin = config%getFloat('n_twin')
177  prm%spr = config%getFloat('s_pr')
178  prm%h0_TwinTwin = config%getFloat('h0_twintwin')
179 
180  ! expand: family => system
181  xi_twin_0 = math_expand(xi_twin_0,n_tw)
182 
183  ! sanity checks
184  if (prm%gdot0_twin <= 0.0_preal) extmsg = trim(extmsg)//' gdot0_twin'
185  if (prm%n_twin <= 0.0_preal) extmsg = trim(extmsg)//' n_twin'
186 
187  else twinactive
188  xi_twin_0 = emptyrealarray
189  allocate(prm%gamma_twin_char,source=emptyrealarray)
190  allocate(prm%interaction_TwinTwin(0,0))
191  endif twinactive
192 
193 !--------------------------------------------------------------------------------------------------
194 ! slip-twin related parameters
195  slipandtwinactive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
196  prm%h0_TwinSlip = config%getFloat('h0_twinslip')
197  prm%interaction_SlipTwin = lattice_interaction_slipbytwin(n_sl,n_tw,&
198  config%getFloats('interaction_sliptwin'), &
199  config%getString('lattice_structure'))
200  prm%interaction_TwinSlip = lattice_interaction_twinbyslip(n_tw,n_sl,&
201  config%getFloats('interaction_twinslip'), &
202  config%getString('lattice_structure'))
203  else slipandtwinactive
204  allocate(prm%interaction_SlipTwin(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0
205  allocate(prm%interaction_TwinSlip(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0
206  prm%h0_TwinSlip = 0.0_preal
207  endif slipandtwinactive
208 
209 !--------------------------------------------------------------------------------------------------
210 ! output pararameters
211  prm%output = config%getStrings('(output)',defaultval=emptystringarray)
212 
213 !--------------------------------------------------------------------------------------------------
214 ! allocate state arrays
215  nipcmyphase = count(material_phaseat == p) * discretization_nip
216  sizedotstate = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
217  + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
218  sizestate = sizedotstate
219 
220  call material_allocateplasticstate(p,nipcmyphase,sizestate,sizedotstate,0)
221 
222 !--------------------------------------------------------------------------------------------------
223 ! state aliases and initialization
224  startindex = 1
225  endindex = prm%sum_N_sl
226  stt%xi_slip => plasticstate(p)%state (startindex:endindex,:)
227  stt%xi_slip = spread(xi_slip_0, 2, nipcmyphase)
228  dot%xi_slip => plasticstate(p)%dotState(startindex:endindex,:)
229  plasticstate(p)%atol(startindex:endindex) = config%getFloat('atol_xi',defaultval=1.0_preal)
230  if(any(plasticstate(p)%atol(startindex:endindex) < 0.0_preal)) extmsg = trim(extmsg)//' atol_xi'
231 
232  startindex = endindex + 1
233  endindex = endindex + prm%sum_N_tw
234  stt%xi_twin => plasticstate(p)%state (startindex:endindex,:)
235  stt%xi_twin = spread(xi_twin_0, 2, nipcmyphase)
236  dot%xi_twin => plasticstate(p)%dotState(startindex:endindex,:)
237  plasticstate(p)%atol(startindex:endindex) = config%getFloat('atol_xi',defaultval=1.0_preal)
238  if(any(plasticstate(p)%atol(startindex:endindex) < 0.0_preal)) extmsg = trim(extmsg)//' atol_xi'
239 
240  startindex = endindex + 1
241  endindex = endindex + prm%sum_N_sl
242  stt%gamma_slip => plasticstate(p)%state (startindex:endindex,:)
243  dot%gamma_slip => plasticstate(p)%dotState(startindex:endindex,:)
244  plasticstate(p)%atol(startindex:endindex) = config%getFloat('atol_gamma',defaultval=1.0e-6_preal)
245  if(any(plasticstate(p)%atol(startindex:endindex) < 0.0_preal)) extmsg = trim(extmsg)//' atol_gamma'
246  ! global alias
247  plasticstate(p)%slipRate => plasticstate(p)%dotState(startindex:endindex,:)
248 
249  startindex = endindex + 1
250  endindex = endindex + prm%sum_N_tw
251  stt%gamma_twin => plasticstate(p)%state (startindex:endindex,:)
252  dot%gamma_twin => plasticstate(p)%dotState(startindex:endindex,:)
253  plasticstate(p)%atol(startindex:endindex) = config%getFloat('atol_gamma',defaultval=1.0e-6_preal)
254  if(any(plasticstate(p)%atol(startindex:endindex) < 0.0_preal)) extmsg = trim(extmsg)//' atol_gamma'
255 
256  plasticstate(p)%state0 = plasticstate(p)%state ! ToDo: this could be done centrally
257 
258  end associate
259 
260 !--------------------------------------------------------------------------------------------------
261 ! exit if any parameter is out of range
262  if (extmsg /= '') call io_error(211,ext_msg=trim(extmsg)//'('//plasticity_phenopowerlaw_label//')')
263 
264  enddo
265 
266 end subroutine plastic_phenopowerlaw_init
267 
268 
269 !--------------------------------------------------------------------------------------------------
272 ! equally (Taylor assumption). Twinning happens only in untwinned volume
273 !--------------------------------------------------------------------------------------------------
274 pure module subroutine plastic_phenopowerlaw_lpanditstangent(lp,dlp_dmp,mp,instance,of)
275 
276  real(pReal), dimension(3,3), intent(out) :: &
277  Lp
278  real(pReal), dimension(3,3,3,3), intent(out) :: &
279  dLp_dMp
280 
281  real(pReal), dimension(3,3), intent(in) :: &
282  Mp
283  integer, intent(in) :: &
284  instance, &
285  of
286 
287  integer :: &
288  i,k,l,m,n
289  real(pReal), dimension(param(instance)%sum_N_sl) :: &
290  gdot_slip_pos,gdot_slip_neg, &
291  dgdot_dtauslip_pos,dgdot_dtauslip_neg
292  real(pReal), dimension(param(instance)%sum_N_tw) :: &
293  gdot_twin,dgdot_dtautwin
294 
295  lp = 0.0_preal
296  dlp_dmp = 0.0_preal
297 
298  associate(prm => param(instance))
299 
300  call kinetics_slip(mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
301  slipsystems: do i = 1, prm%sum_N_sl
302  lp = lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i)
303  forall (k=1:3,l=1:3,m=1:3,n=1:3) &
304  dlp_dmp(k,l,m,n) = dlp_dmp(k,l,m,n) &
305  + dgdot_dtauslip_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) &
306  + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
307  enddo slipsystems
308 
309  call kinetics_twin(mp,instance,of,gdot_twin,dgdot_dtautwin)
310  twinsystems: do i = 1, prm%sum_N_tw
311  lp = lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i)
312  forall (k=1:3,l=1:3,m=1:3,n=1:3) &
313  dlp_dmp(k,l,m,n) = dlp_dmp(k,l,m,n) &
314  + dgdot_dtautwin(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
315  enddo twinsystems
316 
317  end associate
318 
319 end subroutine plastic_phenopowerlaw_lpanditstangent
320 
321 
322 !--------------------------------------------------------------------------------------------------
324 !--------------------------------------------------------------------------------------------------
325 module subroutine plastic_phenopowerlaw_dotstate(mp,instance,of)
326 
327  real(pReal), dimension(3,3), intent(in) :: &
328  Mp
329  integer, intent(in) :: &
330  instance, &
331  of
332 
333  real(pReal) :: &
334  c_SlipSlip,c_TwinSlip,c_TwinTwin, &
335  xi_slip_sat_offset,&
336  sumGamma,sumF
337  real(pReal), dimension(param(instance)%sum_N_sl) :: &
338  left_SlipSlip,right_SlipSlip, &
339  gdot_slip_pos,gdot_slip_neg
340 
341  associate(prm => param(instance), stt => state(instance), dot => dotstate(instance))
342 
343  sumgamma = sum(stt%gamma_slip(:,of))
344  sumf = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)
345 
346 !--------------------------------------------------------------------------------------------------
347 ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
348  c_slipslip = prm%h0_slipslip * (1.0_preal + prm%c_1*sumf** prm%c_2)
349  c_twinslip = prm%h0_TwinSlip * sumgamma**prm%c_3
350  c_twintwin = prm%h0_TwinTwin * sumf**prm%c_4
351 
352 !--------------------------------------------------------------------------------------------------
353 ! calculate left and right vectors
354  left_slipslip = 1.0_preal + prm%H_int
355  xi_slip_sat_offset = prm%spr*sqrt(sumf)
356  right_slipslip = abs(1.0_preal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip &
357  * sign(1.0_preal,1.0_preal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset))
358 
359 !--------------------------------------------------------------------------------------------------
360 ! shear rates
361  call kinetics_slip(mp,instance,of,gdot_slip_pos,gdot_slip_neg)
362  dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg)
363  call kinetics_twin(mp,instance,of,dot%gamma_twin(:,of))
364 
365 !--------------------------------------------------------------------------------------------------
366 ! hardening
367  dot%xi_slip(:,of) = c_slipslip * left_slipslip * &
368  matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_slipslip) &
369  + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of))
370 
371  dot%xi_twin(:,of) = c_twinslip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) &
372  + c_twintwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of))
373  end associate
374 
375 end subroutine plastic_phenopowerlaw_dotstate
376 
377 
378 !--------------------------------------------------------------------------------------------------
380 !--------------------------------------------------------------------------------------------------
381 module subroutine plastic_phenopowerlaw_results(instance,group)
382 
383  integer, intent(in) :: instance
384  character(len=*), intent(in) :: group
385 
386  integer :: o
387 
388  associate(prm => param(instance), stt => state(instance))
389  outputsloop: do o = 1,size(prm%output)
390  select case(trim(prm%output(o)))
391 
392  case('resistance_slip')
393  if(prm%sum_N_sl>0) call results_writedataset(group,stt%xi_slip, 'xi_sl', &
394  'resistance against plastic slip','Pa')
395  case('accumulatedshear_slip')
396  if(prm%sum_N_sl>0) call results_writedataset(group,stt%gamma_slip,'gamma_sl', &
397  'plastic shear','1')
398 
399  case('resistance_twin')
400  if(prm%sum_N_tw>0) call results_writedataset(group,stt%xi_twin, 'xi_tw', &
401  'resistance against twinning','Pa')
402  case('accumulatedshear_twin')
403  if(prm%sum_N_tw>0) call results_writedataset(group,stt%gamma_twin,'gamma_tw', &
404  'twinning shear','1')
405 
406  end select
407  enddo outputsloop
408  end associate
409 
410 end subroutine plastic_phenopowerlaw_results
411 
412 
413 !--------------------------------------------------------------------------------------------------
415 ! stress.
417 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
418 ! have the optional arguments at the end.
419 !--------------------------------------------------------------------------------------------------
420 pure subroutine kinetics_slip(Mp,instance,of, &
421  gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
422 
423  real(pReal), dimension(3,3), intent(in) :: &
424  Mp
425  integer, intent(in) :: &
426  instance, &
427  of
428 
429  real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: &
430  gdot_slip_pos, &
431  gdot_slip_neg
432  real(pReal), intent(out), optional, dimension(param(instance)%sum_N_sl) :: &
433  dgdot_dtau_slip_pos, &
434  dgdot_dtau_slip_neg
435 
436  real(pReal), dimension(param(instance)%sum_N_sl) :: &
437  tau_slip_pos, &
438  tau_slip_neg
439  integer :: i
440 
441  associate(prm => param(instance), stt => state(instance))
442 
443  do i = 1, prm%sum_N_sl
444  tau_slip_pos(i) = math_tensordot(mp,prm%nonSchmid_pos(1:3,1:3,i))
445  tau_slip_neg(i) = merge(math_tensordot(mp,prm%nonSchmid_neg(1:3,1:3,i)), &
446  0.0_preal, prm%nonSchmidActive)
447  enddo
448 
449  where(dneq0(tau_slip_pos))
450  gdot_slip_pos = prm%gdot0_slip * merge(0.5_preal,1.0_preal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
451  * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos)
452  else where
453  gdot_slip_pos = 0.0_preal
454  end where
455 
456  where(dneq0(tau_slip_neg))
457  gdot_slip_neg = prm%gdot0_slip * 0.5_preal & ! only used if non-Schmid active, always 1/2
458  * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg)
459  else where
460  gdot_slip_neg = 0.0_preal
461  end where
462 
463  if (present(dgdot_dtau_slip_pos)) then
464  where(dneq0(gdot_slip_pos))
465  dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos
466  else where
467  dgdot_dtau_slip_pos = 0.0_preal
468  end where
469  endif
470  if (present(dgdot_dtau_slip_neg)) then
471  where(dneq0(gdot_slip_neg))
472  dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg
473  else where
474  dgdot_dtau_slip_neg = 0.0_preal
475  end where
476  endif
477  end associate
478 
479 end subroutine kinetics_slip
480 
481 
482 !--------------------------------------------------------------------------------------------------
484 ! stress. Twinning is assumed to take place only in untwinned volume.
486 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
487 ! have the optional arguments at the end.
488 !--------------------------------------------------------------------------------------------------
489 pure subroutine kinetics_twin(Mp,instance,of,&
490  gdot_twin,dgdot_dtau_twin)
491 
492  real(pReal), dimension(3,3), intent(in) :: &
493  Mp
494  integer, intent(in) :: &
495  instance, &
496  of
497 
498  real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: &
499  gdot_twin
500  real(pReal), dimension(param(instance)%sum_N_tw), intent(out), optional :: &
501  dgdot_dtau_twin
502 
503  real(pReal), dimension(param(instance)%sum_N_tw) :: &
504  tau_twin
505  integer :: i
506 
507  associate(prm => param(instance), stt => state(instance))
508 
509  do i = 1, prm%sum_N_tw
510  tau_twin(i) = math_tensordot(mp,prm%P_tw(1:3,1:3,i))
511  enddo
512 
513  where(tau_twin > 0.0_preal)
514  gdot_twin = (1.0_preal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction
515  * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
516  else where
517  gdot_twin = 0.0_preal
518  end where
519 
520  if (present(dgdot_dtau_twin)) then
521  where(dneq0(gdot_twin))
522  dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin
523  else where
524  dgdot_dtau_twin = 0.0_preal
525  end where
526  endif
527 
528  end associate
529 
530 end subroutine kinetics_twin
531 
532 end submodule plastic_phenopowerlaw
config
Reads in the material configuration from file.
Definition: config.f90:13
constitutive
elasticity, plasticity, internal microstructure state
Definition: constitutive.f90:10