DAMASK with grid solvers  Revision: v2.0.3-2204-gdb1f2151
The Düsseldorf Advanced Material Simulation Kit with Grid Solvers
crystallite.f90
Go to the documentation of this file.
1 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/crystallite.f90"
2 # 1 "<built-in>"
3 # 1 "<command-line>"
4 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/crystallite.f90"
5 !--------------------------------------------------------------------------------------------------
13 !--------------------------------------------------------------------------------------------------
14 
15 module crystallite
16  use prec
17  use io
18  use hdf5_utilities
19  use damask_interface
20  use config
21  use debug
22  use numerics
23  use rotations
24  use math
25  use fesolving
26  use material
27  use constitutive
28  use discretization
29  use lattice
30  use results
31 
32  implicit none
33  private
34 
35  real(preal), dimension(:,:,:), allocatable, public :: &
36  crystallite_dt
37  real(preal), dimension(:,:,:), allocatable :: &
38  crystallite_subdt, & !< substepped time increment of each grain
39  crystallite_subfrac, & !< already calculated fraction of increment
40  crystallite_substep
41  type(rotation), dimension(:,:,:), allocatable :: &
42  crystallite_orientation
43  real(preal), dimension(:,:,:,:,:), allocatable, public, protected :: &
44  crystallite_fe, & !< current "elastic" def grad (end of converged time step)
45  crystallite_p, & !< 1st Piola-Kirchhoff stress per grain
46  crystallite_s0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
47  crystallite_fp0, & !< plastic def grad at start of FE inc
48  crystallite_fi0, & !< intermediate def grad at start of FE inc
49  crystallite_f0, & !< def grad at start of FE inc
50  crystallite_lp0, & !< plastic velocitiy grad at start of FE inc
51  crystallite_li0
52  real(preal), dimension(:,:,:,:,:), allocatable, public :: &
53  crystallite_s, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
54  crystallite_partioneds0, & !< 2nd Piola-Kirchhoff stress vector at start of homog inc
55  crystallite_fp, & !< current plastic def grad (end of converged time step)
56  crystallite_partionedfp0,& !< plastic def grad at start of homog inc
57  crystallite_fi, & !< current intermediate def grad (end of converged time step)
58  crystallite_partionedfi0,& !< intermediate def grad at start of homog inc
59  crystallite_partionedf, & !< def grad to be reached at end of homog inc
60  crystallite_partionedf0, & !< def grad at start of homog inc
61  crystallite_lp, & !< current plastic velocitiy grad (end of converged time step)
62  crystallite_partionedlp0, & !< plastic velocity grad at start of homog inc
63  crystallite_li, & !< current intermediate velocitiy grad (end of converged time step)
64  crystallite_partionedli0
65  real(preal), dimension(:,:,:,:,:), allocatable :: &
66  crystallite_subfp0,& !< plastic def grad at start of crystallite inc
67  crystallite_subfi0,& !< intermediate def grad at start of crystallite inc
68  crystallite_subf, & !< def grad to be reached at end of crystallite inc
69  crystallite_subf0, & !< def grad at start of crystallite inc
70  crystallite_sublp0,& !< plastic velocity grad at start of crystallite inc
71  crystallite_subli0
72  real(preal), dimension(:,:,:,:,:,:,:), allocatable, public, protected :: &
73  crystallite_dpdf
74  logical, dimension(:,:,:), allocatable, public :: &
75  crystallite_requested
76  logical, dimension(:,:,:), allocatable :: &
77  crystallite_converged, & !< convergence flag
78  crystallite_todo, & !< flag to indicate need for further computation
79  crystallite_localplasticity
80 
81  type :: toutput
82  character(len=pStringLen), allocatable, dimension(:) :: &
83  label
84  end type toutput
85  type(toutput), allocatable, dimension(:) :: output_constituent
86 
87  type :: tnumerics
88  integer :: &
89  ijacolpresiduum, & !< frequency of Jacobian update of residuum in Lp
90  nstate, & !< state loop limit
91  nstress
92  real(preal) :: &
93  substepmincryst, & !< minimum (relative) size of sub-step allowed during cutback
94  substepsizecryst, & !< size of first substep when cutback
95  substepsizelp, & !< size of first substep when cutback in Lp calculation
96  substepsizeli, & !< size of first substep when cutback in Li calculation
97  stepincreasecryst, & !< increase of next substep size when previous substep converged
98  rtol_crystallitestate, & !< relative tolerance in state loop
99  rtol_crystallitestress, & !< relative tolerance in stress loop
100  atol_crystallitestress
101  end type tnumerics
102 
103  type(tnumerics) :: num ! numerics parameters. Better name?
104 
105  procedure(), pointer :: integratestate
106 
107  public :: &
108  crystallite_init, &
109  crystallite_stress, &
110  crystallite_stresstangent, &
111  crystallite_orientations, &
112  crystallite_push33toref, &
113  crystallite_results, &
114  crystallite_restartwrite, &
115  crystallite_restartread, &
116  crystallite_forward
117 
118 contains
119 
120 
121 !--------------------------------------------------------------------------------------------------
123 !--------------------------------------------------------------------------------------------------
124 subroutine crystallite_init
125 
126  logical, dimension(discretization_nIP,discretization_nElem) :: devnull
127  integer :: &
128  c, & !< counter in integration point component loop
129  i, & !< counter in integration point loop
130  e, & !< counter in element loop
131  cmax, & !< maximum number of integration point components
132  imax, & !< maximum number of integration points
133  emax, & !< maximum number of elements
134  myncomponents
135 
136  write(6,'(/,a)') ' <<<+- crystallite init -+>>>'
137 
138  cmax = homogenization_maxngrains
139  imax = discretization_nip
140  emax = discretization_nelem
141 
142  allocate(crystallite_partionedf(3,3,cmax,imax,emax),source=0.0_preal)
143 
144  allocate(crystallite_s0, &
145  crystallite_f0, crystallite_fi0,crystallite_fp0, &
146  crystallite_li0,crystallite_lp0, &
147  crystallite_partioneds0, &
148  crystallite_partionedf0,crystallite_partionedfp0,crystallite_partionedfi0, &
149  crystallite_partionedlp0,crystallite_partionedli0, &
150  crystallite_s,crystallite_p, &
151  crystallite_fe,crystallite_fi,crystallite_fp, &
152  crystallite_li,crystallite_lp, &
153  crystallite_subf,crystallite_subf0, &
154  crystallite_subfp0,crystallite_subfi0, &
155  crystallite_subli0,crystallite_sublp0, &
156  source = crystallite_partionedf)
157 
158  allocate(crystallite_dpdf(3,3,3,3,cmax,imax,emax),source=0.0_preal)
159 
160  allocate(crystallite_dt(cmax,imax,emax),source=0.0_preal)
161  allocate(crystallite_subdt,crystallite_subfrac,crystallite_substep, &
162  source = crystallite_dt)
163 
164  allocate(crystallite_orientation(cmax,imax,emax))
165 
166  allocate(crystallite_localplasticity(cmax,imax,emax), source=.true.)
167  allocate(crystallite_requested(cmax,imax,emax), source=.false.)
168  allocate(crystallite_todo(cmax,imax,emax), source=.false.)
169  allocate(crystallite_converged(cmax,imax,emax), source=.true.)
170 
171  num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultval=1.0e-3_preal)
172  num%subStepSizeCryst = config_numerics%getFloat('substepsizecryst', defaultval=0.25_preal)
173  num%stepIncreaseCryst = config_numerics%getFloat('stepincreasecryst', defaultval=1.5_preal)
174 
175  num%subStepSizeLp = config_numerics%getFloat('substepsizelp', defaultval=0.5_preal)
176  num%subStepSizeLi = config_numerics%getFloat('substepsizeli', defaultval=0.5_preal)
177 
178  num%rtol_crystalliteState = config_numerics%getFloat('rtol_crystallitestate', defaultval=1.0e-6_preal)
179  num%rtol_crystalliteStress = config_numerics%getFloat('rtol_crystallitestress',defaultval=1.0e-6_preal)
180  num%atol_crystalliteStress = config_numerics%getFloat('atol_crystallitestress',defaultval=1.0e-8_preal)
181 
182  num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultval=1)
183 
184  num%nState = config_numerics%getInt ('nstate', defaultval=20)
185  num%nStress = config_numerics%getInt ('nstress', defaultval=40)
186 
187  if(num%subStepMinCryst <= 0.0_preal) call io_error(301,ext_msg='subStepMinCryst')
188  if(num%subStepSizeCryst <= 0.0_preal) call io_error(301,ext_msg='subStepSizeCryst')
189  if(num%stepIncreaseCryst <= 0.0_preal) call io_error(301,ext_msg='stepIncreaseCryst')
190 
191  if(num%subStepSizeLp <= 0.0_preal) call io_error(301,ext_msg='subStepSizeLp')
192  if(num%subStepSizeLi <= 0.0_preal) call io_error(301,ext_msg='subStepSizeLi')
193 
194  if(num%rtol_crystalliteState <= 0.0_preal) call io_error(301,ext_msg='rtol_crystalliteState')
195  if(num%rtol_crystalliteStress <= 0.0_preal) call io_error(301,ext_msg='rtol_crystalliteStress')
196  if(num%atol_crystalliteStress <= 0.0_preal) call io_error(301,ext_msg='atol_crystalliteStress')
197 
198  if(num%iJacoLpresiduum < 1) call io_error(301,ext_msg='iJacoLpresiduum')
199 
200  if(num%nState < 1) call io_error(301,ext_msg='nState')
201  if(num%nStress< 1) call io_error(301,ext_msg='nStress')
202 
203  select case(numerics_integrator)
204  case(1)
205  integratestate => integratestatefpi
206  case(2)
207  integratestate => integratestateeuler
208  case(3)
209  integratestate => integratestateadaptiveeuler
210  case(4)
211  integratestate => integratestaterk4
212  case(5)
213  integratestate => integratestaterkck45
214  end select
215 
216  allocate(output_constituent(size(config_phase)))
217  do c = 1, size(config_phase)
218 
219  allocate(output_constituent(c)%label(1))
220  output_constituent(c)%label(1)= 'GfortranBug86277'
221  output_constituent(c)%label = config_phase(c)%getStrings('(output)',defaultval=output_constituent(c)%label )
222  if (output_constituent(c)%label (1) == 'GfortranBug86277') output_constituent(c)%label = [character(len=pStringLen)::]
223 
224 
225 
226  enddo
227 
228  call config_deallocate('material.config/phase')
229 
230 !--------------------------------------------------------------------------------------------------
231 ! initialize
232  !$OMP PARALLEL DO PRIVATE(myNcomponents,i,c)
233  do e = fesolving_execelem(1),fesolving_execelem(2)
234  myncomponents = homogenization_ngrains(material_homogenizationat(e))
235  do i = fesolving_execip(1), fesolving_execip(2); do c = 1, myncomponents
236  crystallite_fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! plastic def gradient reflects init orientation
237  crystallite_fp0(1:3,1:3,c,i,e) = crystallite_fp0(1:3,1:3,c,i,e) &
238  / math_det33(crystallite_fp0(1:3,1:3,c,i,e))**(1.0_preal/3.0_preal)
239  crystallite_fi0(1:3,1:3,c,i,e) = constitutive_initialfi(c,i,e)
240  crystallite_f0(1:3,1:3,c,i,e) = math_i3
241  crystallite_localplasticity(c,i,e) = phase_localplasticity(material_phaseat(c,e))
242  crystallite_fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_fi0(1:3,1:3,c,i,e), &
243  crystallite_fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration
244  crystallite_fp(1:3,1:3,c,i,e) = crystallite_fp0(1:3,1:3,c,i,e)
245  crystallite_fi(1:3,1:3,c,i,e) = crystallite_fi0(1:3,1:3,c,i,e)
246  crystallite_requested(c,i,e) = .true.
247  enddo; enddo
248  enddo
249  !$OMP END PARALLEL DO
250 
251  if(any(.not. crystallite_localplasticity) .and. .not. usepingpong) call io_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal?
252 
253  crystallite_partionedfp0 = crystallite_fp0
254  crystallite_partionedfi0 = crystallite_fi0
255  crystallite_partionedf0 = crystallite_f0
256  crystallite_partionedf = crystallite_f0
257 
258  call crystallite_orientations()
259 
260  !$OMP PARALLEL DO
261  do e = fesolving_execelem(1),fesolving_execelem(2)
262  do i = fesolving_execip(1),fesolving_execip(2)
263  do c = 1,homogenization_ngrains(material_homogenizationat(e))
264  call constitutive_dependentstate(crystallite_partionedf0(1:3,1:3,c,i,e), &
265  crystallite_partionedfp0(1:3,1:3,c,i,e), &
266  c,i,e) ! update dependent state variables to be consistent with basic states
267  enddo
268  enddo
269  enddo
270  !$OMP END PARALLEL DO
271 
272  devnull = crystallite_stress()
273  call crystallite_stresstangent
274 
275 # 283 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/crystallite.f90"
276 
277 end subroutine crystallite_init
278 
279 
280 !--------------------------------------------------------------------------------------------------
282 !--------------------------------------------------------------------------------------------------
283 function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
284 
285  logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress
286  real(preal), intent(in), optional :: &
287  dummyargumenttopreventinternalcompilererrorwithgcc
288  real(preal) :: &
289  formersubstep
290  integer :: &
291  niterationcrystallite, & ! number of iterations in crystallite loop
292  c, & !< counter in integration point component loop
293  i, & !< counter in integration point loop
294  e, & !< counter in element loop
295  startip, endip, &
296  s
297 
298 # 325 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/crystallite.f90"
299 
300 !--------------------------------------------------------------------------------------------------
301 ! initialize to starting condition
302  crystallite_substep = 0.0_preal
303  !$OMP PARALLEL DO
304  elementlooping1: do e = fesolving_execelem(1),fesolving_execelem(2)
305  do i = fesolving_execip(1),fesolving_execip(2); do c = 1,homogenization_ngrains(material_homogenizationat(e))
306  homogenizationrequestscalculation: if (crystallite_requested(c,i,e)) then
307  plasticstate(material_phaseat(c,e))%subState0( :,material_phasememberat(c,i,e)) = &
308  plasticstate(material_phaseat(c,e))%partionedState0(:,material_phasememberat(c,i,e))
309 
310  do s = 1, phase_nsources(material_phaseat(c,e))
311  sourcestate(material_phaseat(c,e))%p(s)%subState0( :,material_phasememberat(c,i,e)) = &
312  sourcestate(material_phaseat(c,e))%p(s)%partionedState0(:,material_phasememberat(c,i,e))
313  enddo
314  crystallite_subfp0(1:3,1:3,c,i,e) = crystallite_partionedfp0(1:3,1:3,c,i,e)
315  crystallite_sublp0(1:3,1:3,c,i,e) = crystallite_partionedlp0(1:3,1:3,c,i,e)
316  crystallite_subfi0(1:3,1:3,c,i,e) = crystallite_partionedfi0(1:3,1:3,c,i,e)
317  crystallite_subli0(1:3,1:3,c,i,e) = crystallite_partionedli0(1:3,1:3,c,i,e)
318  crystallite_subf0(1:3,1:3,c,i,e) = crystallite_partionedf0(1:3,1:3,c,i,e)
319  crystallite_subfrac(c,i,e) = 0.0_preal
320  crystallite_substep(c,i,e) = 1.0_preal/num%subStepSizeCryst
321  crystallite_todo(c,i,e) = .true.
322  crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst
323  endif homogenizationrequestscalculation
324  enddo; enddo
325  enddo elementlooping1
326  !$OMP END PARALLEL DO
327 
328  singlerun: if (fesolving_execelem(1) == fesolving_execelem(2) .and. &
329  fesolving_execip(1) == fesolving_execip(2)) then
330  startip = fesolving_execip(1)
331  endip = startip
332  else singlerun
333  startip = 1
334  endip = discretization_nip
335  endif singlerun
336 
337  niterationcrystallite = 0
338  cutbacklooping: do while (any(crystallite_todo(:,startip:endip,fesolving_execelem(1):fesolving_execelem(2))))
339  niterationcrystallite = niterationcrystallite + 1
340 
341 
342 
343 
344 
345  !$OMP PARALLEL DO PRIVATE(formerSubStep)
346  elementlooping3: do e = fesolving_execelem(1),fesolving_execelem(2)
347  do i = fesolving_execip(1),fesolving_execip(2)
348  do c = 1,homogenization_ngrains(material_homogenizationat(e))
349 !--------------------------------------------------------------------------------------------------
350 ! wind forward
351  if (crystallite_converged(c,i,e)) then
352  formersubstep = crystallite_substep(c,i,e)
353  crystallite_subfrac(c,i,e) = crystallite_subfrac(c,i,e) + crystallite_substep(c,i,e)
354  crystallite_substep(c,i,e) = min(1.0_preal - crystallite_subfrac(c,i,e), &
355  num%stepIncreaseCryst * crystallite_substep(c,i,e))
356 
357  crystallite_todo(c,i,e) = crystallite_substep(c,i,e) > 0.0_preal ! still time left to integrate on?
358  if (crystallite_todo(c,i,e)) then
359  crystallite_subf0(1:3,1:3,c,i,e) = crystallite_subf(1:3,1:3,c,i,e)
360  crystallite_sublp0(1:3,1:3,c,i,e) = crystallite_lp(1:3,1:3,c,i,e)
361  crystallite_subli0(1:3,1:3,c,i,e) = crystallite_li(1:3,1:3,c,i,e)
362  crystallite_subfp0(1:3,1:3,c,i,e) = crystallite_fp(1:3,1:3,c,i,e)
363  crystallite_subfi0(1:3,1:3,c,i,e) = crystallite_fi(1:3,1:3,c,i,e)
364  !if abbrevation, make c and p private in omp
365  plasticstate( material_phaseat(c,e))%subState0(:,material_phasememberat(c,i,e)) &
366  = plasticstate(material_phaseat(c,e))%state( :,material_phasememberat(c,i,e))
367  do s = 1, phase_nsources(material_phaseat(c,e))
368  sourcestate( material_phaseat(c,e))%p(s)%subState0(:,material_phasememberat(c,i,e)) &
369  = sourcestate(material_phaseat(c,e))%p(s)%state( :,material_phasememberat(c,i,e))
370  enddo
371  endif
372 
373 !--------------------------------------------------------------------------------------------------
374 ! cut back (reduced time and restore)
375  else
376  crystallite_substep(c,i,e) = num%subStepSizeCryst * crystallite_substep(c,i,e)
377  crystallite_fp(1:3,1:3,c,i,e) = crystallite_subfp0(1:3,1:3,c,i,e)
378  crystallite_fi(1:3,1:3,c,i,e) = crystallite_subfi0(1:3,1:3,c,i,e)
379  crystallite_s(1:3,1:3,c,i,e) = crystallite_s0(1:3,1:3,c,i,e)
380  if (crystallite_substep(c,i,e) < 1.0_preal) then ! actual (not initial) cutback
381  crystallite_lp(1:3,1:3,c,i,e) = crystallite_sublp0(1:3,1:3,c,i,e)
382  crystallite_li(1:3,1:3,c,i,e) = crystallite_subli0(1:3,1:3,c,i,e)
383  endif
384  plasticstate(material_phaseat(c,e))%state( :,material_phasememberat(c,i,e)) &
385  = plasticstate(material_phaseat(c,e))%subState0(:,material_phasememberat(c,i,e))
386  do s = 1, phase_nsources(material_phaseat(c,e))
387  sourcestate( material_phaseat(c,e))%p(s)%state( :,material_phasememberat(c,i,e)) &
388  = sourcestate(material_phaseat(c,e))%p(s)%subState0(:,material_phasememberat(c,i,e))
389  enddo
390 
391  ! cant restore dotState here, since not yet calculated in first cutback after initialization
392  crystallite_todo(c,i,e) = crystallite_substep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair)
393  endif
394 
395 !--------------------------------------------------------------------------------------------------
396 ! prepare for integration
397  if (crystallite_todo(c,i,e)) then
398  crystallite_subf(1:3,1:3,c,i,e) = crystallite_subf0(1:3,1:3,c,i,e) &
399  + crystallite_substep(c,i,e) *( crystallite_partionedf(1:3,1:3,c,i,e) &
400  -crystallite_partionedf0(1:3,1:3,c,i,e))
401  crystallite_fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subf(1:3,1:3,c,i,e), &
402  math_inv33(crystallite_fp(1:3,1:3,c,i,e))), &
403  math_inv33(crystallite_fi(1:3,1:3,c,i,e)))
404  crystallite_subdt(c,i,e) = crystallite_substep(c,i,e) * crystallite_dt(c,i,e)
405  crystallite_converged(c,i,e) = .false.
406  endif
407 
408  enddo
409  enddo
410  enddo elementlooping3
411  !$OMP END PARALLEL DO
412 
413 !--------------------------------------------------------------------------------------------------
414 ! integrate --- requires fully defined state array (basic + dependent state)
415  if (any(crystallite_todo)) call integratestate ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation
416  where(.not. crystallite_converged .and. crystallite_substep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further
417  crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
418 
419 
420  enddo cutbacklooping
421 
422 ! return whether converged or not
423  crystallite_stress = .false.
424  elementlooping5: do e = fesolving_execelem(1),fesolving_execelem(2)
425  do i = fesolving_execip(1),fesolving_execip(2)
426  crystallite_stress(i,e) = all(crystallite_converged(:,i,e))
427  enddo
428  enddo elementlooping5
429 
430 end function crystallite_stress
431 
432 
433 !--------------------------------------------------------------------------------------------------
435 !--------------------------------------------------------------------------------------------------
436 subroutine crystallite_stresstangent
437 
438  integer :: &
439  c, & !< counter in integration point component loop
440  i, & !< counter in integration point loop
441  e, & !< counter in element loop
442  o, &
443  p
444 
445  real(preal), dimension(3,3) :: devnull, &
446  invsubfp0,invsubfi0,invfp,invfi, &
447  temp_33_1, temp_33_2, temp_33_3, temp_33_4
448  real(preal), dimension(3,3,3,3) :: dsdfe, &
449  dsdf, &
450  dsdfi, &
451  dlids, & ! tangent in lattice configuration
452  dlidfi, &
453  dlpds, &
454  dlpdfi, &
455  dfids, &
456  dfpinvdf, &
457  rhs_3333, &
458  lhs_3333, &
459  temp_3333
460  real(preal), dimension(9,9):: temp_99
461  logical :: error
462 
463  !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,o,p, &
464  !$OMP invSubFp0,invSubFi0,invFp,invFi, &
465  !$OMP rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error)
466  elementlooping: do e = fesolving_execelem(1),fesolving_execelem(2)
467  do i = fesolving_execip(1),fesolving_execip(2)
468  do c = 1,homogenization_ngrains(material_homogenizationat(e))
469 
470  call constitutive_sanditstangents(devnull,dsdfe,dsdfi, &
471  crystallite_fe(1:3,1:3,c,i,e), &
472  crystallite_fi(1:3,1:3,c,i,e),c,i,e)
473  call constitutive_lianditstangents(devnull,dlids,dlidfi, &
474  crystallite_s(1:3,1:3,c,i,e), &
475  crystallite_fi(1:3,1:3,c,i,e), &
476  c,i,e)
477 
478  invfp = math_inv33(crystallite_fp(1:3,1:3,c,i,e))
479  invfi = math_inv33(crystallite_fi(1:3,1:3,c,i,e))
480  invsubfp0 = math_inv33(crystallite_subfp0(1:3,1:3,c,i,e))
481  invsubfi0 = math_inv33(crystallite_subfi0(1:3,1:3,c,i,e))
482 
483  if (sum(abs(dlids)) < tol_math_check) then
484  dfids = 0.0_preal
485  else
486  lhs_3333 = 0.0_preal; rhs_3333 = 0.0_preal
487  do o=1,3; do p=1,3
488  lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
489  + crystallite_subdt(c,i,e)*matmul(invsubfi0,dlidfi(1:3,1:3,o,p))
490  lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
491  + invfi*invfi(p,o)
492  rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
493  - crystallite_subdt(c,i,e)*matmul(invsubfi0,dlids(1:3,1:3,o,p))
494  enddo; enddo
495  call math_invert(temp_99,error,math_3333to99(lhs_3333))
496  if (error) then
497  call io_warning(warning_id=600,el=e,ip=i,g=c, &
498  ext_msg='inversion error in analytic tangent calculation')
499  dfids = 0.0_preal
500  else
501  dfids = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
502  endif
503  dlids = math_mul3333xx3333(dlidfi,dfids) + dlids
504  endif
505 
506  call constitutive_lpanditstangents(devnull,dlpds,dlpdfi, &
507  crystallite_s(1:3,1:3,c,i,e), &
508  crystallite_fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration
509  dlpds = math_mul3333xx3333(dlpdfi,dfids) + dlpds
510 
511 !--------------------------------------------------------------------------------------------------
512 ! calculate dSdF
513  temp_33_1 = transpose(matmul(invfp,invfi))
514  temp_33_2 = matmul(crystallite_subf(1:3,1:3,c,i,e),invsubfp0)
515  temp_33_3 = matmul(matmul(crystallite_subf(1:3,1:3,c,i,e),invfp), invsubfi0)
516 
517  do o=1,3; do p=1,3
518  rhs_3333(p,o,1:3,1:3) = matmul(dsdfe(p,o,1:3,1:3),temp_33_1)
519  temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dlpds(1:3,1:3,p,o)), invfi) &
520  + matmul(temp_33_3,dlids(1:3,1:3,p,o))
521  enddo; enddo
522  lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dsdfe,temp_3333) &
523  + math_mul3333xx3333(dsdfi,dfids)
524 
525  call math_invert(temp_99,error,math_identity2nd(9)+math_3333to99(lhs_3333))
526  if (error) then
527  call io_warning(warning_id=600,el=e,ip=i,g=c, &
528  ext_msg='inversion error in analytic tangent calculation')
529  dsdf = rhs_3333
530  else
531  dsdf = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
532  endif
533 
534 !--------------------------------------------------------------------------------------------------
535 ! calculate dFpinvdF
536  temp_3333 = math_mul3333xx3333(dlpds,dsdf)
537  do o=1,3; do p=1,3
538  dfpinvdf(1:3,1:3,p,o) = -crystallite_subdt(c,i,e) &
539  * matmul(invsubfp0, matmul(temp_3333(1:3,1:3,p,o),invfi))
540  enddo; enddo
541 
542 !--------------------------------------------------------------------------------------------------
543 ! assemble dPdF
544  temp_33_1 = matmul(crystallite_s(1:3,1:3,c,i,e),transpose(invfp))
545  temp_33_2 = matmul(invfp,temp_33_1)
546  temp_33_3 = matmul(crystallite_subf(1:3,1:3,c,i,e),invfp)
547  temp_33_4 = matmul(temp_33_3,crystallite_s(1:3,1:3,c,i,e))
548 
549  crystallite_dpdf(1:3,1:3,1:3,1:3,c,i,e) = 0.0_preal
550  do p=1,3
551  crystallite_dpdf(p,1:3,p,1:3,c,i,e) = transpose(temp_33_2)
552  enddo
553  do o=1,3; do p=1,3
554  crystallite_dpdf(1:3,1:3,p,o,c,i,e) = crystallite_dpdf(1:3,1:3,p,o,c,i,e) &
555  + matmul(matmul(crystallite_subf(1:3,1:3,c,i,e), &
556  dfpinvdf(1:3,1:3,p,o)),temp_33_1) &
557  + matmul(matmul(temp_33_3,dsdf(1:3,1:3,p,o)), &
558  transpose(invfp)) &
559  + matmul(temp_33_4,transpose(dfpinvdf(1:3,1:3,p,o)))
560  enddo; enddo
561 
562  enddo; enddo
563  enddo elementlooping
564  !$OMP END PARALLEL DO
565 
566 end subroutine crystallite_stresstangent
567 
568 
569 !--------------------------------------------------------------------------------------------------
571 !--------------------------------------------------------------------------------------------------
572 subroutine crystallite_orientations
573 
574  integer &
575  c, & !< counter in integration point component loop
576  i, & !< counter in integration point loop
577  e
578 
579  !$OMP PARALLEL DO
580  do e = fesolving_execelem(1),fesolving_execelem(2)
581  do i = fesolving_execip(1),fesolving_execip(2)
582  do c = 1,homogenization_ngrains(material_homogenizationat(e))
583  call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalpart(crystallite_fe(1:3,1:3,c,i,e))))
584  enddo; enddo; enddo
585  !$OMP END PARALLEL DO
586 
587  nonlocalpresent: if (any(plasticstate%nonLocal)) then
588  !$OMP PARALLEL DO
589  do e = fesolving_execelem(1),fesolving_execelem(2)
590  do i = fesolving_execip(1),fesolving_execip(2)
591  if (plasticstate(material_phaseat(1,e))%nonLocal) &
592  call plastic_nonlocal_updatecompatibility(crystallite_orientation, &
593  phase_plasticityinstance(material_phaseat(i,e)),i,e)
594  enddo; enddo
595  !$OMP END PARALLEL DO
596  endif nonlocalpresent
597 
598 end subroutine crystallite_orientations
599 
600 
601 !--------------------------------------------------------------------------------------------------
603 !--------------------------------------------------------------------------------------------------
604 function crystallite_push33toref(ipc,ip,el, tensor33)
605 
606  real(preal), dimension(3,3) :: crystallite_push33toref
607  real(preal), dimension(3,3), intent(in) :: tensor33
608  real(preal), dimension(3,3) :: t
609  integer, intent(in):: &
610  el, &
611  ip, &
612  ipc
613 
614  t = matmul(material_orientation0(ipc,ip,el)%asMatrix(), & ! ToDo: initial orientation correct?
615  transpose(math_inv33(crystallite_subf(1:3,1:3,ipc,ip,el))))
616  crystallite_push33toref = matmul(transpose(t),matmul(tensor33,t))
617 
618 end function crystallite_push33toref
619 
620 
621 !--------------------------------------------------------------------------------------------------
623 !--------------------------------------------------------------------------------------------------
624 subroutine crystallite_results
625 
626  integer :: p,o
627  real(preal), allocatable, dimension(:,:,:) :: selected_tensors
628  type(rotation), allocatable, dimension(:) :: selected_rotations
629  character(len=pStringLen) :: group,structurelabel
630 
631  do p=1,size(config_name_phase)
632  group = trim('current/constituent')//'/'//trim(config_name_phase(p))//'/generic'
633 
634  call results_closegroup(results_addgroup(group))
635 
636  do o = 1, size(output_constituent(p)%label)
637  select case (output_constituent(p)%label(o))
638  case('f')
639  selected_tensors = select_tensors(crystallite_partionedf,p)
640  call results_writedataset(group,selected_tensors,'F',&
641  'deformation gradient','1')
642  case('fe')
643  selected_tensors = select_tensors(crystallite_fe,p)
644  call results_writedataset(group,selected_tensors,'Fe',&
645  'elastic deformation gradient','1')
646  case('fp')
647  selected_tensors = select_tensors(crystallite_fp,p)
648  call results_writedataset(group,selected_tensors,'Fp',&
649  'plastic deformation gradient','1')
650  case('fi')
651  selected_tensors = select_tensors(crystallite_fi,p)
652  call results_writedataset(group,selected_tensors,'Fi',&
653  'inelastic deformation gradient','1')
654  case('lp')
655  selected_tensors = select_tensors(crystallite_lp,p)
656  call results_writedataset(group,selected_tensors,'Lp',&
657  'plastic velocity gradient','1/s')
658  case('li')
659  selected_tensors = select_tensors(crystallite_li,p)
660  call results_writedataset(group,selected_tensors,'Li',&
661  'inelastic velocity gradient','1/s')
662  case('p')
663  selected_tensors = select_tensors(crystallite_p,p)
664  call results_writedataset(group,selected_tensors,'P',&
665  'First Piola-Kirchoff stress','Pa')
666  case('s')
667  selected_tensors = select_tensors(crystallite_s,p)
668  call results_writedataset(group,selected_tensors,'S',&
669  'Second Piola-Kirchoff stress','Pa')
670  case('orientation')
671  select case(lattice_structure(p))
672  case(lattice_iso_id)
673  structurelabel = 'iso'
674  case(lattice_fcc_id)
675  structurelabel = 'fcc'
676  case(lattice_bcc_id)
677  structurelabel = 'bcc'
678  case(lattice_bct_id)
679  structurelabel = 'bct'
680  case(lattice_hex_id)
681  structurelabel = 'hex'
682  case(lattice_ort_id)
683  structurelabel = 'ort'
684  end select
685  selected_rotations = select_rotations(crystallite_orientation,p)
686  call results_writedataset(group,selected_rotations,'orientation',&
687  'crystal orientation as quaternion',structurelabel)
688  end select
689  enddo
690  enddo
691 
692  contains
693 
694  !------------------------------------------------------------------------------------------------
696  !------------------------------------------------------------------------------------------------
697  function select_tensors(dataset,instance)
698 
699  integer, intent(in) :: instance
700  real(preal), dimension(:,:,:,:,:), intent(in) :: dataset
701  real(preal), allocatable, dimension(:,:,:) :: select_tensors
702  integer :: e,i,c,j
703 
704  allocate(select_tensors(3,3,count(material_phaseat==instance)*discretization_nip))
705 
706  j=0
707  do e = 1, size(material_phaseat,2)
708  do i = 1, discretization_nip
709  do c = 1, size(material_phaseat,1) !ToDo: this needs to be changed for varying Ngrains
710  if (material_phaseat(c,e) == instance) then
711  j = j + 1
712  select_tensors(1:3,1:3,j) = dataset(1:3,1:3,c,i,e)
713  endif
714  enddo
715  enddo
716  enddo
717 
718  end function select_tensors
719 
720 
721 !--------------------------------------------------------------------------------------------------
723 !--------------------------------------------------------------------------------------------------
724  function select_rotations(dataset,instance)
725 
726  integer, intent(in) :: instance
727  type(rotation), dimension(:,:,:), intent(in) :: dataset
728  type(rotation), allocatable, dimension(:) :: select_rotations
729  integer :: e,i,c,j
730 
731  allocate(select_rotations(count(material_phaseat==instance)*homogenization_maxngrains*discretization_nip))
732 
733  j=0
734  do e = 1, size(material_phaseat,2)
735  do i = 1, discretization_nip
736  do c = 1, size(material_phaseat,1) !ToDo: this needs to be changed for varying Ngrains
737  if (material_phaseat(c,e) == instance) then
738  j = j + 1
739  select_rotations(j) = dataset(c,i,e)
740  endif
741  enddo
742  enddo
743  enddo
744 
745  end function select_rotations
746 
747 end subroutine crystallite_results
748 
749 
750 !--------------------------------------------------------------------------------------------------
753 !--------------------------------------------------------------------------------------------------
754 logical function integratestress(ipc,ip,el,timeFraction)
755 
756  integer, intent(in):: el, & ! element index
757  ip, & ! integration point index
758  ipc ! grain index
759  real(preal), optional, intent(in) :: timefraction ! fraction of timestep
760 
761  real(preal), dimension(3,3):: f, & ! deformation gradient at end of timestep
762  fp_new, & ! plastic deformation gradient at end of timestep
763  invfp_new, & ! inverse of Fp_new
764  invfp_current, & ! inverse of Fp_current
765  lpguess, & ! current guess for plastic velocity gradient
766  lpguess_old, & ! known last good guess for plastic velocity gradient
767  lp_constitutive, & ! plastic velocity gradient resulting from constitutive law
768  residuumlp, & ! current residuum of plastic velocity gradient
769  residuumlp_old, & ! last residuum of plastic velocity gradient
770  deltalp, & ! direction of next guess
771  fi_new, & ! gradient of intermediate deformation stages
772  invfi_new, &
773  invfi_current, & ! inverse of Fi_current
774  liguess, & ! current guess for intermediate velocity gradient
775  liguess_old, & ! known last good guess for intermediate velocity gradient
776  li_constitutive, & ! intermediate velocity gradient resulting from constitutive law
777  residuumli, & ! current residuum of intermediate velocity gradient
778  residuumli_old, & ! last residuum of intermediate velocity gradient
779  deltali, & ! direction of next guess
780  fe, & ! elastic deformation gradient
781  fe_new, &
782  s, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
783  a, &
784  b, &
785  temp_33
786  real(preal), dimension(9) :: temp_9 ! needed for matrix inversion by LAPACK
787  integer, dimension(9) :: devnull_9 ! needed for matrix inversion by LAPACK
788  real(preal), dimension(9,9) :: drlp_dlp, & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme)
789  drli_dli ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme)
790  real(preal), dimension(3,3,3,3):: ds_dfe, & ! partial derivative of 2nd Piola-Kirchhoff stress
791  ds_dfi, &
792  dfe_dlp, & ! partial derivative of elastic deformation gradient
793  dfe_dli, &
794  dfi_dli, &
795  dlp_dfi, &
796  dli_dfi, &
797  dlp_ds, &
798  dli_ds
799  real(preal) steplengthlp, &
800  steplengthli, &
801  dt, & ! time increment
802  atol_lp, &
803  atol_li, &
804  devnull
805  integer niterationstresslp, & ! number of stress integrations
806  niterationstressli, & ! number of inner stress integrations
807  ierr, & ! error indicator for LAPACK
808  o, &
809  p, &
810  jacocounterlp, &
811  jacocounterli ! counters to check for Jacobian update
812  logical :: error
813  external :: &
814  dgesv
815 
816  integratestress = .false.
817 
818  if (present(timefraction)) then
819  dt = crystallite_subdt(ipc,ip,el) * timefraction
820  f = crystallite_subf0(1:3,1:3,ipc,ip,el) &
821  + (crystallite_subf(1:3,1:3,ipc,ip,el) - crystallite_subf0(1:3,1:3,ipc,ip,el)) * timefraction
822  else
823  dt = crystallite_subdt(ipc,ip,el)
824  f = crystallite_subf(1:3,1:3,ipc,ip,el)
825  endif
826 
827  lpguess = crystallite_lp(1:3,1:3,ipc,ip,el) ! take as first guess
828  liguess = crystallite_li(1:3,1:3,ipc,ip,el) ! take as first guess
829 
830  call math_invert33(invfp_current,devnull,error,crystallite_subfp0(1:3,1:3,ipc,ip,el))
831  if (error) return ! error
832  call math_invert33(invfi_current,devnull,error,crystallite_subfi0(1:3,1:3,ipc,ip,el))
833  if (error) return ! error
834 
835  a = matmul(f,invfp_current) ! intermediate tensor needed later to calculate dFe_dLp
836 
837  jacocounterli = 0
838  steplengthli = 1.0_preal
839  residuumli_old = 0.0_preal
840  liguess_old = liguess
841 
842  niterationstressli = 0
843  liloop: do
844  niterationstressli = niterationstressli + 1
845  if (niterationstressli>num%nStress) return ! error
846 
847  invfi_new = matmul(invfi_current,math_i3 - dt*liguess)
848  fi_new = math_inv33(invfi_new)
849 
850  jacocounterlp = 0
851  steplengthlp = 1.0_preal
852  residuumlp_old = 0.0_preal
853  lpguess_old = lpguess
854 
855  niterationstresslp = 0
856  lploop: do
857  niterationstresslp = niterationstresslp + 1
858  if (niterationstresslp>num%nStress) return ! error
859 
860  b = math_i3 - dt*lpguess
861  fe = matmul(matmul(a,b), invfi_new)
862  call constitutive_sanditstangents(s, ds_dfe, ds_dfi, &
863  fe, fi_new, ipc, ip, el)
864 
865  call constitutive_lpanditstangents(lp_constitutive, dlp_ds, dlp_dfi, &
866  s, fi_new, ipc, ip, el)
867 
868  !* update current residuum and check for convergence of loop
869  atol_lp = max(num%rtol_crystalliteStress * max(norm2(lpguess),norm2(lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
870  num%atol_crystalliteStress) ! minimum lower cutoff
871  residuumlp = lpguess - lp_constitutive
872 
873  if (any(ieee_is_nan(residuumlp))) then
874  return ! error
875  elseif (norm2(residuumlp) < atol_lp) then ! converged if below absolute tolerance
876  exit lploop
877  elseif (niterationstresslp == 1 .or. norm2(residuumlp) < norm2(residuumlp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
878  residuumlp_old = residuumlp ! ...remember old values and...
879  lpguess_old = lpguess
880  steplengthlp = 1.0_preal ! ...proceed with normal step length (calculate new search direction)
881  else ! not converged and residuum not improved...
882  steplengthlp = num%subStepSizeLp * steplengthlp ! ...try with smaller step length in same direction
883  lpguess = lpguess_old &
884  + deltalp * steplengthlp
885  cycle lploop
886  endif
887 
888  !* calculate Jacobian for correction term
889  if (mod(jacocounterlp, num%iJacoLpresiduum) == 0) then
890  jacocounterlp = jacocounterlp + 1
891 
892  do o=1,3; do p=1,3
893  dfe_dlp(o,1:3,p,1:3) = a(o,p)*transpose(invfi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
894  enddo; enddo
895  dfe_dlp = - dt * dfe_dlp
896  drlp_dlp = math_identity2nd(9) &
897  - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dlp_ds,ds_dfe),dfe_dlp))
898  temp_9 = math_33to9(residuumlp)
899  call dgesv(9,1,drlp_dlp,9,devnull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
900  if (ierr /= 0) return ! error
901  deltalp = - math_9to33(temp_9)
902  endif
903 
904  lpguess = lpguess &
905  + deltalp * steplengthlp
906  enddo lploop
907 
908  call constitutive_lianditstangents(li_constitutive, dli_ds, dli_dfi, &
909  s, fi_new, ipc, ip, el)
910 
911  !* update current residuum and check for convergence of loop
912  atol_li = max(num%rtol_crystalliteStress * max(norm2(liguess),norm2(li_constitutive)), & ! absolute tolerance from largest acceptable relative error
913  num%atol_crystalliteStress) ! minimum lower cutoff
914  residuumli = liguess - li_constitutive
915  if (any(ieee_is_nan(residuumli))) then
916  return ! error
917  elseif (norm2(residuumli) < atol_li) then ! converged if below absolute tolerance
918  exit liloop
919  elseif (niterationstressli == 1 .or. norm2(residuumli) < norm2(residuumli_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
920  residuumli_old = residuumli ! ...remember old values and...
921  liguess_old = liguess
922  steplengthli = 1.0_preal ! ...proceed with normal step length (calculate new search direction)
923  else ! not converged and residuum not improved...
924  steplengthli = num%subStepSizeLi * steplengthli ! ...try with smaller step length in same direction
925  liguess = liguess_old &
926  + deltali * steplengthli
927  cycle liloop
928  endif
929 
930  !* calculate Jacobian for correction term
931  if (mod(jacocounterli, num%iJacoLpresiduum) == 0) then
932  jacocounterli = jacocounterli + 1
933 
934  temp_33 = matmul(matmul(a,b),invfi_current)
935  do o=1,3; do p=1,3
936  dfe_dli(1:3,o,1:3,p) = -dt*math_i3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
937  dfi_dli(1:3,o,1:3,p) = -dt*math_i3(o,p)*invfi_current
938  enddo; enddo
939  do o=1,3; do p=1,3
940  dfi_dli(1:3,1:3,o,p) = matmul(matmul(fi_new,dfi_dli(1:3,1:3,o,p)),fi_new)
941  enddo; enddo
942  drli_dli = math_identity2nd(9) &
943  - math_3333to99(math_mul3333xx3333(dli_ds, math_mul3333xx3333(ds_dfe, dfe_dli) &
944  + math_mul3333xx3333(ds_dfi, dfi_dli))) &
945  - math_3333to99(math_mul3333xx3333(dli_dfi, dfi_dli))
946  temp_9 = math_33to9(residuumli)
947  call dgesv(9,1,drli_dli,9,devnull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
948  if (ierr /= 0) return ! error
949  deltali = - math_9to33(temp_9)
950  endif
951 
952  liguess = liguess &
953  + deltali * steplengthli
954  enddo liloop
955 
956  invfp_new = matmul(invfp_current,b)
957  call math_invert33(fp_new,devnull,error,invfp_new)
958  if (error) return ! error
959  fp_new = fp_new / math_det33(fp_new)**(1.0_preal/3.0_preal) ! regularize
960  fe_new = matmul(matmul(f,invfp_new),invfi_new)
961 
962  integratestress = .true.
963  crystallite_p(1:3,1:3,ipc,ip,el) = matmul(matmul(f,invfp_new),matmul(s,transpose(invfp_new)))
964  crystallite_s(1:3,1:3,ipc,ip,el) = s
965  crystallite_lp(1:3,1:3,ipc,ip,el) = lpguess
966  crystallite_li(1:3,1:3,ipc,ip,el) = liguess
967  crystallite_fp(1:3,1:3,ipc,ip,el) = fp_new
968  crystallite_fi(1:3,1:3,ipc,ip,el) = fi_new
969  crystallite_fe(1:3,1:3,ipc,ip,el) = fe_new
970 
971 end function integratestress
972 
973 
974 !--------------------------------------------------------------------------------------------------
977 !--------------------------------------------------------------------------------------------------
978 subroutine integratestatefpi
979 
980  integer :: &
981  NiterationState, & !< number of iterations in state loop
982  e, & !< element index in element loop
983  i, & !< integration point index in ip loop
984  g, & !< grain index in grain loop
985  p, &
986  c, &
987  s, &
988  sizeDotState
989  real(pReal) :: &
990  zeta
991  real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: &
992  residuum_plastic ! residuum for plastic state
993  real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
994  residuum_source ! residuum for source state
995  logical :: &
996  nonlocalBroken
997 
998  nonlocalbroken = .false.
999  !$OMP PARALLEL DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c)
1000  do e = fesolving_execelem(1),fesolving_execelem(2)
1001  do i = fesolving_execip(1),fesolving_execip(2)
1002  do g = 1,homogenization_ngrains(material_homogenizationat(e))
1003  if(crystallite_todo(g,i,e) .and. (.not. nonlocalbroken .or. crystallite_localplasticity(g,i,e)) ) then
1004 
1005  p = material_phaseat(g,e); c = material_phasememberat(g,i,e)
1006 
1007  call constitutive_collectdotstate(crystallite_s(1:3,1:3,g,i,e), &
1008  crystallite_partionedf0, &
1009  crystallite_fi(1:3,1:3,g,i,e), &
1010  crystallite_partionedfp0, &
1011  crystallite_subdt(g,i,e), g,i,e)
1012  crystallite_todo(g,i,e) = all(.not. ieee_is_nan(plasticstate(p)%dotState(:,c)))
1013  do s = 1, phase_nsources(p)
1014  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. ieee_is_nan(sourcestate(p)%p(s)%dotState(:,c)))
1015  enddo
1016  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1017  nonlocalbroken = .true.
1018  if(.not. crystallite_todo(g,i,e)) cycle
1019 
1020  sizedotstate = plasticstate(p)%sizeDotState
1021  plasticstate(p)%state(1:sizedotstate,c) = plasticstate(p)%subState0(1:sizedotstate,c) &
1022  + plasticstate(p)%dotState (1:sizedotstate,c) &
1023  * crystallite_subdt(g,i,e)
1024  do s = 1, phase_nsources(p)
1025  sizedotstate = sourcestate(p)%p(s)%sizeDotState
1026  sourcestate(p)%p(s)%state(1:sizedotstate,c) = sourcestate(p)%p(s)%subState0(1:sizedotstate,c) &
1027  + sourcestate(p)%p(s)%dotState (1:sizedotstate,c) &
1028  * crystallite_subdt(g,i,e)
1029  enddo
1030 
1031  iteration: do niterationstate = 1, num%nState
1032 
1033  plasticstate(p)%previousDotState2(:,c) = merge(plasticstate(p)%previousDotState(:,c),&
1034  0.0_preal,&
1035  niterationstate > 1)
1036  plasticstate(p)%previousDotState (:,c) = plasticstate(p)%dotState(:,c)
1037  do s = 1, phase_nsources(p)
1038  sourcestate(p)%p(s)%previousDotState2(:,c) = merge(sourcestate(p)%p(s)%previousDotState(:,c),&
1039  0.0_preal, &
1040  niterationstate > 1)
1041  sourcestate(p)%p(s)%previousDotState (:,c) = sourcestate(p)%p(s)%dotState(:,c)
1042  enddo
1043 
1044  call constitutive_dependentstate(crystallite_partionedf(1:3,1:3,g,i,e), &
1045  crystallite_fp(1:3,1:3,g,i,e), &
1046  g, i, e)
1047 
1048  crystallite_todo(g,i,e) = integratestress(g,i,e)
1049  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1050  nonlocalbroken = .true.
1051  if(.not. crystallite_todo(g,i,e)) exit iteration
1052 
1053  call constitutive_collectdotstate(crystallite_s(1:3,1:3,g,i,e), &
1054  crystallite_partionedf0, &
1055  crystallite_fi(1:3,1:3,g,i,e), &
1056  crystallite_partionedfp0, &
1057  crystallite_subdt(g,i,e), g,i,e)
1058  crystallite_todo(g,i,e) = all(.not. ieee_is_nan(plasticstate(p)%dotState(:,c)))
1059  do s = 1, phase_nsources(p)
1060  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. ieee_is_nan(sourcestate(p)%p(s)%dotState(:,c)))
1061  enddo
1062  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1063  nonlocalbroken = .true.
1064  if(.not. crystallite_todo(g,i,e)) exit iteration
1065 
1066  sizedotstate = plasticstate(p)%sizeDotState
1067  zeta = damper(plasticstate(p)%dotState (:,c), &
1068  plasticstate(p)%previousDotState (:,c), &
1069  plasticstate(p)%previousDotState2(:,c))
1070  plasticstate(p)%dotState(:,c) = plasticstate(p)%dotState(:,c) * zeta &
1071  + plasticstate(p)%previousDotState(:,c) * (1.0_preal - zeta)
1072  residuum_plastic(1:sizedotstate) = plasticstate(p)%state (1:sizedotstate,c) &
1073  - plasticstate(p)%subState0(1:sizedotstate,c) &
1074  - plasticstate(p)%dotState (1:sizedotstate,c) &
1075  * crystallite_subdt(g,i,e)
1076  plasticstate(p)%state(1:sizedotstate,c) = plasticstate(p)%state(1:sizedotstate,c) &
1077  - residuum_plastic(1:sizedotstate)
1078  crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizedotstate), &
1079  plasticstate(p)%state(1:sizedotstate,c), &
1080  plasticstate(p)%atol(1:sizedotstate))
1081  do s = 1, phase_nsources(p)
1082  sizedotstate = sourcestate(p)%p(s)%sizeDotState
1083  zeta = damper(sourcestate(p)%p(s)%dotState (:,c), &
1084  sourcestate(p)%p(s)%previousDotState (:,c), &
1085  sourcestate(p)%p(s)%previousDotState2(:,c))
1086  sourcestate(p)%p(s)%dotState(:,c) = sourcestate(p)%p(s)%dotState(:,c) * zeta &
1087  + sourcestate(p)%p(s)%previousDotState(:,c)* (1.0_preal - zeta)
1088  residuum_source(1:sizedotstate) = sourcestate(p)%p(s)%state (1:sizedotstate,c) &
1089  - sourcestate(p)%p(s)%subState0(1:sizedotstate,c) &
1090  - sourcestate(p)%p(s)%dotState (1:sizedotstate,c) &
1091  * crystallite_subdt(g,i,e)
1092  sourcestate(p)%p(s)%state(1:sizedotstate,c) = sourcestate(p)%p(s)%state(1:sizedotstate,c) &
1093  - residuum_source(1:sizedotstate)
1094  crystallite_converged(g,i,e) = &
1095  crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizedotstate), &
1096  sourcestate(p)%p(s)%state(1:sizedotstate,c), &
1097  sourcestate(p)%p(s)%atol(1:sizedotstate))
1098  enddo
1099 
1100  if(crystallite_converged(g,i,e)) then
1101  crystallite_todo(g,i,e) = statejump(g,i,e)
1102  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1103  nonlocalbroken = .true.
1104  exit iteration
1105  endif
1106 
1107  enddo iteration
1108 
1109  endif
1110  enddo; enddo; enddo
1111  !$OMP END PARALLEL DO
1112 
1113  if(nonlocalbroken) where(.not. crystallite_localplasticity) crystallite_todo = .false.
1114  if (any(plasticstate(:)%nonlocal)) call nonlocalconvergencecheck
1115 
1116  contains
1117 
1118  !--------------------------------------------------------------------------------------------------
1120  !--------------------------------------------------------------------------------------------------
1121  real(pReal) pure function damper(current,previous,previous2)
1122 
1123  real(preal), dimension(:), intent(in) ::&
1124  current, previous, previous2
1125 
1126  real(preal) :: dot_prod12, dot_prod22
1127 
1128  dot_prod12 = dot_product(current - previous, previous - previous2)
1129  dot_prod22 = dot_product(previous - previous2, previous - previous2)
1130  if ((dot_product(current,previous) < 0.0_preal .or. dot_prod12 < 0.0_preal) .and. dot_prod22 > 0.0_preal) then
1131  damper = 0.75_preal + 0.25_preal * tanh(2.0_preal + 4.0_preal * dot_prod12 / dot_prod22)
1132  else
1133  damper = 1.0_preal
1134  endif
1135 
1136  end function damper
1137 
1138 end subroutine integratestatefpi
1139 
1140 
1141 !--------------------------------------------------------------------------------------------------
1143 !--------------------------------------------------------------------------------------------------
1144 subroutine integratestateeuler
1145 
1146  integer :: &
1147  e, & !< element index in element loop
1148  i, & !< integration point index in ip loop
1149  g, & !< grain index in grain loop
1150  p, &
1151  c, &
1152  s, &
1153  sizeDotState
1154  logical :: &
1155  nonlocalBroken
1156 
1157  nonlocalbroken = .false.
1158  !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c)
1159  do e = fesolving_execelem(1),fesolving_execelem(2)
1160  do i = fesolving_execip(1),fesolving_execip(2)
1161  do g = 1,homogenization_ngrains(material_homogenizationat(e))
1162  if(crystallite_todo(g,i,e) .and. (.not. nonlocalbroken .or. crystallite_localplasticity(g,i,e)) ) then
1163 
1164  p = material_phaseat(g,e); c = material_phasememberat(g,i,e)
1165 
1166  call constitutive_collectdotstate(crystallite_s(1:3,1:3,g,i,e), &
1167  crystallite_partionedf0, &
1168  crystallite_fi(1:3,1:3,g,i,e), &
1169  crystallite_partionedfp0, &
1170  crystallite_subdt(g,i,e), g,i,e)
1171  crystallite_todo(g,i,e) = all(.not. ieee_is_nan(plasticstate(p)%dotState(:,c)))
1172  do s = 1, phase_nsources(p)
1173  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. ieee_is_nan(sourcestate(p)%p(s)%dotState(:,c)))
1174  enddo
1175  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1176  nonlocalbroken = .true.
1177  if(.not. crystallite_todo(g,i,e)) cycle
1178 
1179  sizedotstate = plasticstate(p)%sizeDotState
1180  plasticstate(p)%state(1:sizedotstate,c) = plasticstate(p)%subState0(1:sizedotstate,c) &
1181  + plasticstate(p)%dotState (1:sizedotstate,c) &
1182  * crystallite_subdt(g,i,e)
1183  do s = 1, phase_nsources(p)
1184  sizedotstate = sourcestate(p)%p(s)%sizeDotState
1185  sourcestate(p)%p(s)%state(1:sizedotstate,c) = sourcestate(p)%p(s)%subState0(1:sizedotstate,c) &
1186  + sourcestate(p)%p(s)%dotState (1:sizedotstate,c) &
1187  * crystallite_subdt(g,i,e)
1188  enddo
1189 
1190  crystallite_todo(g,i,e) = statejump(g,i,e)
1191  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1192  nonlocalbroken = .true.
1193  if(.not. crystallite_todo(g,i,e)) cycle
1194 
1195  call constitutive_dependentstate(crystallite_partionedf(1:3,1:3,g,i,e), &
1196  crystallite_fp(1:3,1:3,g,i,e), &
1197  g, i, e)
1198 
1199  crystallite_todo(g,i,e) = integratestress(g,i,e)
1200  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1201  nonlocalbroken = .true.
1202 
1203  crystallite_converged(g,i,e) = crystallite_todo(g,i,e)
1204 
1205  endif
1206  enddo; enddo; enddo
1207  !$OMP END PARALLEL DO
1208 
1209  if(nonlocalbroken) where(.not. crystallite_localplasticity) crystallite_todo = .false.
1210  if (any(plasticstate(:)%nonlocal)) call nonlocalconvergencecheck
1211 
1212 end subroutine integratestateeuler
1213 
1214 
1215 !--------------------------------------------------------------------------------------------------
1217 !--------------------------------------------------------------------------------------------------
1218 subroutine integratestateadaptiveeuler
1219 
1220  integer :: &
1221  e, & ! element index in element loop
1222  i, & ! integration point index in ip loop
1223  g, & ! grain index in grain loop
1224  p, &
1225  c, &
1226  s, &
1227  sizeDotState
1228  logical :: &
1229  nonlocalBroken
1230 
1231  real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
1232  residuum_plastic
1233  real(pReal), dimension(constitutive_source_maxSizeDotState,& maxval(phase_Nsources), & homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
1234  residuum_source
1235 
1236 
1237  nonlocalbroken = .false.
1238  !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
1239  do e = fesolving_execelem(1),fesolving_execelem(2)
1240  do i = fesolving_execip(1),fesolving_execip(2)
1241  do g = 1,homogenization_ngrains(material_homogenizationat(e))
1242  if(crystallite_todo(g,i,e) .and. (.not. nonlocalbroken .or. crystallite_localplasticity(g,i,e)) ) then
1243 
1244  p = material_phaseat(g,e); c = material_phasememberat(g,i,e)
1245 
1246  call constitutive_collectdotstate(crystallite_s(1:3,1:3,g,i,e), &
1247  crystallite_partionedf0, &
1248  crystallite_fi(1:3,1:3,g,i,e), &
1249  crystallite_partionedfp0, &
1250  crystallite_subdt(g,i,e), g,i,e)
1251  crystallite_todo(g,i,e) = all(.not. ieee_is_nan(plasticstate(p)%dotState(:,c)))
1252  do s = 1, phase_nsources(p)
1253  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. ieee_is_nan(sourcestate(p)%p(s)%dotState(:,c)))
1254  enddo
1255  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1256  nonlocalbroken = .true.
1257  if(.not. crystallite_todo(g,i,e)) cycle
1258 
1259  sizedotstate = plasticstate(p)%sizeDotState
1260 
1261  residuum_plastic(1:sizedotstate,g,i,e) = plasticstate(p)%dotstate(1:sizedotstate,c) &
1262  * (- 0.5_preal * crystallite_subdt(g,i,e))
1263  plasticstate(p)%state(1:sizedotstate,c) = plasticstate(p)%subState0(1:sizedotstate,c) &
1264  + plasticstate(p)%dotstate(1:sizedotstate,c) * crystallite_subdt(g,i,e)
1265  do s = 1, phase_nsources(p)
1266  sizedotstate = sourcestate(p)%p(s)%sizeDotState
1267 
1268  residuum_source(1:sizedotstate,s,g,i,e) = sourcestate(p)%p(s)%dotstate(1:sizedotstate,c) &
1269  * (- 0.5_preal * crystallite_subdt(g,i,e))
1270  sourcestate(p)%p(s)%state(1:sizedotstate,c) = sourcestate(p)%p(s)%subState0(1:sizedotstate,c) &
1271  + sourcestate(p)%p(s)%dotstate(1:sizedotstate,c) * crystallite_subdt(g,i,e)
1272  enddo
1273 
1274  crystallite_todo(g,i,e) = statejump(g,i,e)
1275  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1276  nonlocalbroken = .true.
1277  if(.not. crystallite_todo(g,i,e)) cycle
1278 
1279  call constitutive_dependentstate(crystallite_partionedf(1:3,1:3,g,i,e), &
1280  crystallite_fp(1:3,1:3,g,i,e), &
1281  g, i, e)
1282 
1283  crystallite_todo(g,i,e) = integratestress(g,i,e)
1284  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1285  nonlocalbroken = .true.
1286  if(.not. crystallite_todo(g,i,e)) cycle
1287 
1288  call constitutive_collectdotstate(crystallite_s(1:3,1:3,g,i,e), &
1289  crystallite_partionedf0, &
1290  crystallite_fi(1:3,1:3,g,i,e), &
1291  crystallite_partionedfp0, &
1292  crystallite_subdt(g,i,e), g,i,e)
1293  crystallite_todo(g,i,e) = all(.not. ieee_is_nan(plasticstate(p)%dotState(:,c)))
1294  do s = 1, phase_nsources(p)
1295  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. ieee_is_nan(sourcestate(p)%p(s)%dotState(:,c)))
1296  enddo
1297  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1298  nonlocalbroken = .true.
1299  if(.not. crystallite_todo(g,i,e)) cycle
1300 
1301 
1302  sizedotstate = plasticstate(p)%sizeDotState
1303 
1304  residuum_plastic(1:sizedotstate,g,i,e) = residuum_plastic(1:sizedotstate,g,i,e) &
1305  + 0.5_preal * plasticstate(p)%dotState(:,c) * crystallite_subdt(g,i,e)
1306 
1307  crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizedotstate,g,i,e), &
1308  plasticstate(p)%state(1:sizedotstate,c), &
1309  plasticstate(p)%atol(1:sizedotstate))
1310 
1311  do s = 1, phase_nsources(p)
1312  sizedotstate = sourcestate(p)%p(s)%sizeDotState
1313 
1314  residuum_source(1:sizedotstate,s,g,i,e) = &
1315  residuum_source(1:sizedotstate,s,g,i,e) + 0.5_preal * sourcestate(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e)
1316 
1317  crystallite_converged(g,i,e) = &
1318  crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizedotstate,s,g,i,e), &
1319  sourcestate(p)%p(s)%state(1:sizedotstate,c), &
1320  sourcestate(p)%p(s)%atol(1:sizedotstate))
1321  enddo
1322 
1323  endif
1324  enddo; enddo; enddo
1325  !$OMP END PARALLEL DO
1326 
1327  if (any(plasticstate(:)%nonlocal)) call nonlocalconvergencecheck
1328 
1329 end subroutine integratestateadaptiveeuler
1330 
1331 
1332 !--------------------------------------------------------------------------------------------------
1334 !--------------------------------------------------------------------------------------------------
1335 subroutine integratestaterk4
1336 
1337  real(pReal), dimension(3,3), parameter :: &
1338  A = reshape([&
1339  0.5_preal, 0.0_preal, 0.0_preal, &
1340  0.0_preal, 0.5_preal, 0.0_preal, &
1341  0.0_preal, 0.0_preal, 1.0_preal], &
1342  [3,3])
1343  real(pReal), dimension(3), parameter :: &
1344  CC = [0.5_preal, 0.5_preal, 1.0_preal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
1345  real(pReal), dimension(4), parameter :: &
1346  B = [1.0_preal/6.0_preal, 1.0_preal/3.0_preal, 1.0_preal/3.0_preal, 1.0_preal/6.0_preal] ! weight of slope used for Runge Kutta integration (final weight divided by 6)
1347 
1348  integer :: &
1349  e, & ! element index in element loop
1350  i, & ! integration point index in ip loop
1351  g, & ! grain index in grain loop
1352  stage, & ! stage index in integration stage loop
1353  n, &
1354  p, &
1355  c, &
1356  s, &
1357  sizeDotState
1358  logical :: &
1359  nonlocalBroken
1360 
1361  nonlocalbroken = .false.
1362  !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
1363  do e = fesolving_execelem(1),fesolving_execelem(2)
1364  do i = fesolving_execip(1),fesolving_execip(2)
1365  do g = 1,homogenization_ngrains(material_homogenizationat(e))
1366  if(crystallite_todo(g,i,e) .and. (.not. nonlocalbroken .or. crystallite_localplasticity(g,i,e)) ) then
1367 
1368  p = material_phaseat(g,e); c = material_phasememberat(g,i,e)
1369 
1370  call constitutive_collectdotstate(crystallite_s(1:3,1:3,g,i,e), &
1371  crystallite_partionedf0, &
1372  crystallite_fi(1:3,1:3,g,i,e), &
1373  crystallite_partionedfp0, &
1374  crystallite_subdt(g,i,e), g,i,e)
1375  crystallite_todo(g,i,e) = all(.not. ieee_is_nan(plasticstate(p)%dotState(:,c)))
1376  do s = 1, phase_nsources(p)
1377  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. ieee_is_nan(sourcestate(p)%p(s)%dotState(:,c)))
1378  enddo
1379  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1380  nonlocalbroken = .true.
1381  if(.not. crystallite_todo(g,i,e)) cycle
1382 
1383  do stage = 1,3
1384 
1385  plasticstate(p)%RK4dotState(stage,:,c) = plasticstate(p)%dotState(:,c)
1386  plasticstate(p)%dotState(:,c) = a(1,stage) * plasticstate(p)%RK4dotState(1,:,c)
1387  do s = 1, phase_nsources(p)
1388  sourcestate(p)%p(s)%RK4dotState(stage,:,c) = sourcestate(p)%p(s)%dotState(:,c)
1389  sourcestate(p)%p(s)%dotState(:,c) = a(1,stage) * sourcestate(p)%p(s)%RK4dotState(1,:,c)
1390  enddo
1391 
1392  do n = 2, stage
1393  plasticstate(p)%dotState(:,c) = plasticstate(p)%dotState(:,c) &
1394  + a(n,stage) * plasticstate(p)%RK4dotState(n,:,c)
1395  do s = 1, phase_nsources(p)
1396  sourcestate(p)%p(s)%dotState(:,c) = sourcestate(p)%p(s)%dotState(:,c) &
1397  + a(n,stage) * sourcestate(p)%p(s)%RK4dotState(n,:,c)
1398  enddo
1399  enddo
1400 
1401  sizedotstate = plasticstate(p)%sizeDotState
1402  plasticstate(p)%state(1:sizedotstate,c) = plasticstate(p)%subState0(1:sizedotstate,c) &
1403  + plasticstate(p)%dotState (1:sizedotstate,c) &
1404  * crystallite_subdt(g,i,e)
1405  do s = 1, phase_nsources(p)
1406  sizedotstate = sourcestate(p)%p(s)%sizeDotState
1407  sourcestate(p)%p(s)%state(1:sizedotstate,c) = sourcestate(p)%p(s)%subState0(1:sizedotstate,c) &
1408  + sourcestate(p)%p(s)%dotState (1:sizedotstate,c) &
1409  * crystallite_subdt(g,i,e)
1410  enddo
1411 
1412  call constitutive_dependentstate(crystallite_partionedf(1:3,1:3,g,i,e), &
1413  crystallite_fp(1:3,1:3,g,i,e), &
1414  g, i, e)
1415 
1416  crystallite_todo(g,i,e) = integratestress(g,i,e,cc(stage))
1417  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1418  nonlocalbroken = .true.
1419  if(.not. crystallite_todo(g,i,e)) exit
1420 
1421  call constitutive_collectdotstate(crystallite_s(1:3,1:3,g,i,e), &
1422  crystallite_partionedf0, &
1423  crystallite_fi(1:3,1:3,g,i,e), &
1424  crystallite_partionedfp0, &
1425  crystallite_subdt(g,i,e)*cc(stage), g,i,e)
1426  crystallite_todo(g,i,e) = all(.not. ieee_is_nan(plasticstate(p)%dotState(:,c)))
1427  do s = 1, phase_nsources(p)
1428  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. ieee_is_nan(sourcestate(p)%p(s)%dotState(:,c)))
1429  enddo
1430  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1431  nonlocalbroken = .true.
1432  if(.not. crystallite_todo(g,i,e)) exit
1433 
1434  enddo
1435 
1436  if(.not. crystallite_todo(g,i,e)) cycle
1437 
1438  sizedotstate = plasticstate(p)%sizeDotState
1439 
1440  plasticstate(p)%RK4dotState(4,:,c) = plasticstate(p)%dotState(:,c)
1441 
1442  plasticstate(p)%dotState(:,c) = matmul(b,plasticstate(p)%RK4dotState(1:4,1:sizedotstate,c))
1443  plasticstate(p)%state(1:sizedotstate,c) = plasticstate(p)%subState0(1:sizedotstate,c) &
1444  + plasticstate(p)%dotState (1:sizedotstate,c) &
1445  * crystallite_subdt(g,i,e)
1446 
1447  do s = 1, phase_nsources(p)
1448  sizedotstate = sourcestate(p)%p(s)%sizeDotState
1449 
1450  sourcestate(p)%p(s)%RK4dotState(4,:,c) = sourcestate(p)%p(s)%dotState(:,c)
1451 
1452  sourcestate(p)%p(s)%dotState(:,c) = matmul(b,sourcestate(p)%p(s)%RK4dotState(1:4,1:sizedotstate,c))
1453  sourcestate(p)%p(s)%state(1:sizedotstate,c) = sourcestate(p)%p(s)%subState0(1:sizedotstate,c) &
1454  + sourcestate(p)%p(s)%dotState (1:sizedotstate,c) &
1455  * crystallite_subdt(g,i,e)
1456  enddo
1457 
1458  crystallite_todo(g,i,e) = statejump(g,i,e)
1459  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1460  nonlocalbroken = .true.
1461  if(.not. crystallite_todo(g,i,e)) cycle
1462 
1463  call constitutive_dependentstate(crystallite_partionedf(1:3,1:3,g,i,e), &
1464  crystallite_fp(1:3,1:3,g,i,e), &
1465  g, i, e)
1466 
1467  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1468  nonlocalbroken = .true.
1469  if(.not. crystallite_todo(g,i,e)) cycle
1470 
1471  crystallite_todo(g,i,e) = integratestress(g,i,e)
1472  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1473  nonlocalbroken = .true.
1474  crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken
1475 
1476  endif
1477  enddo; enddo; enddo
1478  !$OMP END PARALLEL DO
1479 
1480  if(nonlocalbroken) where(.not. crystallite_localplasticity) crystallite_todo = .false.
1481  if (any(plasticstate(:)%nonlocal)) call nonlocalconvergencecheck
1482 
1483 end subroutine integratestaterk4
1484 
1485 
1486 !--------------------------------------------------------------------------------------------------
1489 !--------------------------------------------------------------------------------------------------
1490 subroutine integratestaterkck45
1491 
1492  real(pReal), dimension(5,5), parameter :: &
1493  A = reshape([&
1494  .2_preal, .075_preal, .3_preal, -11.0_preal/54.0_preal, 1631.0_preal/55296.0_preal, &
1495  .0_preal, .225_preal, -.9_preal, 2.5_preal, 175.0_preal/512.0_preal, &
1496  .0_preal, .0_preal, 1.2_preal, -70.0_preal/27.0_preal, 575.0_preal/13824.0_preal, &
1497  .0_preal, .0_preal, .0_preal, 35.0_preal/27.0_preal, 44275.0_preal/110592.0_preal, &
1498  .0_preal, .0_preal, .0_preal, .0_preal, 253.0_preal/4096.0_preal], &
1499  [5,5], order=[2,1])
1500 
1501  real(pReal), dimension(6), parameter :: &
1502  B = &
1503  [37.0_preal/378.0_preal, .0_preal, 250.0_preal/621.0_preal, &
1504  125.0_preal/594.0_preal, .0_preal, 512.0_preal/1771.0_preal], &
1505  db = b - &
1506  [2825.0_preal/27648.0_preal, .0_preal, 18575.0_preal/48384.0_preal,&
1507  13525.0_preal/55296.0_preal, 277.0_preal/14336.0_preal, 0.25_preal]
1508 
1509  real(pReal), dimension(5), parameter :: &
1510  CC = [0.2_preal, 0.3_preal, 0.6_preal, 1.0_preal, 0.875_preal]
1511 
1512  integer :: &
1513  e, & ! element index in element loop
1514  i, & ! integration point index in ip loop
1515  g, & ! grain index in grain loop
1516  stage, & ! stage index in integration stage loop
1517  n, &
1518  p, &
1519  c, &
1520  s, &
1521  sizeDotState
1522  logical :: &
1523  nonlocalBroken
1524 
1525  real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
1526  residuum_plastic
1527  real(pReal), dimension(constitutive_source_maxSizeDotState, & maxval(phase_Nsources), & homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
1528  residuum_source
1529 
1530 
1531  nonlocalbroken = .false.
1532  !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
1533  do e = fesolving_execelem(1),fesolving_execelem(2)
1534  do i = fesolving_execip(1),fesolving_execip(2)
1535  do g = 1,homogenization_ngrains(material_homogenizationat(e))
1536  if(crystallite_todo(g,i,e) .and. (.not. nonlocalbroken .or. crystallite_localplasticity(g,i,e)) ) then
1537 
1538  p = material_phaseat(g,e); c = material_phasememberat(g,i,e)
1539 
1540  call constitutive_collectdotstate(crystallite_s(1:3,1:3,g,i,e), &
1541  crystallite_partionedf0, &
1542  crystallite_fi(1:3,1:3,g,i,e), &
1543  crystallite_partionedfp0, &
1544  crystallite_subdt(g,i,e), g,i,e)
1545  crystallite_todo(g,i,e) = all(.not. ieee_is_nan(plasticstate(p)%dotState(:,c)))
1546  do s = 1, phase_nsources(p)
1547  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. ieee_is_nan(sourcestate(p)%p(s)%dotState(:,c)))
1548  enddo
1549  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1550  nonlocalbroken = .true.
1551  if(.not. crystallite_todo(g,i,e)) cycle
1552 
1553  do stage = 1,5
1554 
1555  plasticstate(p)%RKCK45dotState(stage,:,c) = plasticstate(p)%dotState(:,c)
1556  plasticstate(p)%dotState(:,c) = a(1,stage) * plasticstate(p)%RKCK45dotState(1,:,c)
1557  do s = 1, phase_nsources(p)
1558  sourcestate(p)%p(s)%RKCK45dotState(stage,:,c) = sourcestate(p)%p(s)%dotState(:,c)
1559  sourcestate(p)%p(s)%dotState(:,c) = a(1,stage) * sourcestate(p)%p(s)%RKCK45dotState(1,:,c)
1560  enddo
1561 
1562  do n = 2, stage
1563  plasticstate(p)%dotState(:,c) = plasticstate(p)%dotState(:,c) &
1564  + a(n,stage) * plasticstate(p)%RKCK45dotState(n,:,c)
1565  do s = 1, phase_nsources(p)
1566  sourcestate(p)%p(s)%dotState(:,c) = sourcestate(p)%p(s)%dotState(:,c) &
1567  + a(n,stage) * sourcestate(p)%p(s)%RKCK45dotState(n,:,c)
1568  enddo
1569  enddo
1570 
1571  sizedotstate = plasticstate(p)%sizeDotState
1572  plasticstate(p)%state(1:sizedotstate,c) = plasticstate(p)%subState0(1:sizedotstate,c) &
1573  + plasticstate(p)%dotState (1:sizedotstate,c) &
1574  * crystallite_subdt(g,i,e)
1575  do s = 1, phase_nsources(p)
1576  sizedotstate = sourcestate(p)%p(s)%sizeDotState
1577  sourcestate(p)%p(s)%state(1:sizedotstate,c) = sourcestate(p)%p(s)%subState0(1:sizedotstate,c) &
1578  + sourcestate(p)%p(s)%dotState (1:sizedotstate,c) &
1579  * crystallite_subdt(g,i,e)
1580  enddo
1581 
1582  call constitutive_dependentstate(crystallite_partionedf(1:3,1:3,g,i,e), &
1583  crystallite_fp(1:3,1:3,g,i,e), &
1584  g, i, e)
1585 
1586  crystallite_todo(g,i,e) = integratestress(g,i,e,cc(stage))
1587  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1588  nonlocalbroken = .true.
1589  if(.not. crystallite_todo(g,i,e)) exit
1590 
1591  call constitutive_collectdotstate(crystallite_s(1:3,1:3,g,i,e), &
1592  crystallite_partionedf0, &
1593  crystallite_fi(1:3,1:3,g,i,e), &
1594  crystallite_partionedfp0, &
1595  crystallite_subdt(g,i,e)*cc(stage), g,i,e)
1596  crystallite_todo(g,i,e) = all(.not. ieee_is_nan(plasticstate(p)%dotState(:,c)))
1597  do s = 1, phase_nsources(p)
1598  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. ieee_is_nan(sourcestate(p)%p(s)%dotState(:,c)))
1599  enddo
1600  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1601  nonlocalbroken = .true.
1602  if(.not. crystallite_todo(g,i,e)) exit
1603 
1604  enddo
1605 
1606  if(.not. crystallite_todo(g,i,e)) cycle
1607 
1608  sizedotstate = plasticstate(p)%sizeDotState
1609 
1610  plasticstate(p)%RKCK45dotState(6,:,c) = plasticstate(p)%dotState(:,c)
1611  residuum_plastic(1:sizedotstate,g,i,e) = matmul(db,plasticstate(p)%RKCK45dotState(1:6,1:sizedotstate,c)) &
1612  * crystallite_subdt(g,i,e)
1613  plasticstate(p)%dotState(:,c) = matmul(b,plasticstate(p)%RKCK45dotState(1:6,1:sizedotstate,c))
1614  plasticstate(p)%state(1:sizedotstate,c) = plasticstate(p)%subState0(1:sizedotstate,c) &
1615  + plasticstate(p)%dotState (1:sizedotstate,c) &
1616  * crystallite_subdt(g,i,e)
1617  crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizedotstate,g,i,e), &
1618  plasticstate(p)%state(1:sizedotstate,c), &
1619  plasticstate(p)%atol(1:sizedotstate))
1620 
1621  do s = 1, phase_nsources(p)
1622  sizedotstate = sourcestate(p)%p(s)%sizeDotState
1623 
1624  sourcestate(p)%p(s)%RKCK45dotState(6,:,c) = sourcestate(p)%p(s)%dotState(:,c)
1625  residuum_source(1:sizedotstate,s,g,i,e) = matmul(db,sourcestate(p)%p(s)%RKCK45dotState(1:6,1:sizedotstate,c)) &
1626  * crystallite_subdt(g,i,e)
1627  sourcestate(p)%p(s)%dotState(:,c) = matmul(b,sourcestate(p)%p(s)%RKCK45dotState(1:6,1:sizedotstate,c))
1628  sourcestate(p)%p(s)%state(1:sizedotstate,c) = sourcestate(p)%p(s)%subState0(1:sizedotstate,c) &
1629  + sourcestate(p)%p(s)%dotState (1:sizedotstate,c) &
1630  * crystallite_subdt(g,i,e)
1631  crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. &
1632  converged(residuum_source(1:sizedotstate,s,g,i,e), &
1633  sourcestate(p)%p(s)%state(1:sizedotstate,c), &
1634  sourcestate(p)%p(s)%atol(1:sizedotstate))
1635  enddo
1636  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1637  nonlocalbroken = .true.
1638  if(.not. crystallite_todo(g,i,e)) cycle
1639 
1640  crystallite_todo(g,i,e) = statejump(g,i,e)
1641  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1642  nonlocalbroken = .true.
1643  if(.not. crystallite_todo(g,i,e)) cycle
1644 
1645  call constitutive_dependentstate(crystallite_partionedf(1:3,1:3,g,i,e), &
1646  crystallite_fp(1:3,1:3,g,i,e), &
1647  g, i, e)
1648 
1649  crystallite_todo(g,i,e) = integratestress(g,i,e)
1650  if(.not. (crystallite_todo(g,i,e) .or. crystallite_localplasticity(g,i,e))) &
1651  nonlocalbroken = .true.
1652  crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken
1653 
1654  endif
1655  enddo; enddo; enddo
1656  !$OMP END PARALLEL DO
1657 
1658  if(nonlocalbroken) where(.not. crystallite_localplasticity) crystallite_todo = .false.
1659  if (any(plasticstate(:)%nonlocal)) call nonlocalconvergencecheck
1660 
1661 end subroutine integratestaterkck45
1662 
1663 
1664 !--------------------------------------------------------------------------------------------------
1667 !--------------------------------------------------------------------------------------------------
1668 subroutine nonlocalconvergencecheck
1669 
1670  if (any(.not. crystallite_converged .and. .not. crystallite_localplasticity)) & ! any non-local not yet converged (or broken)...
1671  where( .not. crystallite_localplasticity) crystallite_converged = .false.
1672 
1673 end subroutine nonlocalconvergencecheck
1674 
1675 
1676 !--------------------------------------------------------------------------------------------------
1678 !--------------------------------------------------------------------------------------------------
1679 logical pure function converged(residuum,state,atol)
1680 
1681  real(pReal), intent(in), dimension(:) ::&
1682  residuum, state, atol
1683  real(pReal) :: &
1684  rTol
1685 
1686  rtol = num%rTol_crystalliteState
1687 
1688  converged = all(abs(residuum) <= max(atol, rtol*abs(state)))
1689 
1690 end function converged
1691 
1692 
1693 !--------------------------------------------------------------------------------------------------
1696 !--------------------------------------------------------------------------------------------------
1697 logical function statejump(ipc,ip,el)
1698 
1699  integer, intent(in):: &
1700  el, & ! element index
1701  ip, & ! integration point index
1702  ipc ! grain index
1703 
1704  integer :: &
1705  c, &
1706  p, &
1707  mysource, &
1708  myoffset, &
1709  mysize
1710 
1711  c = material_phasememberat(ipc,ip,el)
1712  p = material_phaseat(ipc,el)
1713 
1714  call constitutive_collectdeltastate(crystallite_s(1:3,1:3,ipc,ip,el), &
1715  crystallite_fe(1:3,1:3,ipc,ip,el), &
1716  crystallite_fi(1:3,1:3,ipc,ip,el), &
1717  ipc,ip,el)
1718 
1719  myoffset = plasticstate(p)%offsetDeltaState
1720  mysize = plasticstate(p)%sizeDeltaState
1721 
1722  if( any(ieee_is_nan(plasticstate(p)%deltaState(1:mysize,c)))) then
1723  statejump = .false.
1724  return
1725  endif
1726 
1727  plasticstate(p)%state(myoffset + 1:myoffset + mysize,c) = &
1728  plasticstate(p)%state(myoffset + 1:myoffset + mysize,c) + plasticstate(p)%deltaState(1:mysize,c)
1729 
1730  do mysource = 1, phase_nsources(p)
1731  myoffset = sourcestate(p)%p(mysource)%offsetDeltaState
1732  mysize = sourcestate(p)%p(mysource)%sizeDeltaState
1733  if (any(ieee_is_nan(sourcestate(p)%p(mysource)%deltaState(1:mysize,c)))) then
1734  statejump = .false.
1735  return
1736  endif
1737  sourcestate(p)%p(mysource)%state(myoffset + 1: myoffset + mysize,c) = &
1738  sourcestate(p)%p(mysource)%state(myoffset + 1: myoffset + mysize,c) + sourcestate(p)%p(mysource)%deltaState(1:mysize,c)
1739  enddo
1740 
1741  statejump = .true.
1742 
1743 end function statejump
1744 
1745 
1746 !--------------------------------------------------------------------------------------------------
1748 ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
1749 !--------------------------------------------------------------------------------------------------
1750 subroutine crystallite_restartwrite
1751 
1752  integer :: i
1753  integer(HID_T) :: filehandle, grouphandle
1754  character(len=pStringLen) :: filename, datasetname
1755 
1756  write(6,'(a)') ' writing field and constitutive data required for restart to file';flush(6)
1757 
1758  write(filename,'(a,i0,a)') trim(getsolverjobname())//'_',worldrank,'.hdf5'
1759  filehandle = hdf5_openfile(filename,'a')
1760 
1761  call hdf5_write(filehandle,crystallite_partionedf,'F')
1762  call hdf5_write(filehandle,crystallite_fp, 'Fp')
1763  call hdf5_write(filehandle,crystallite_fi, 'Fi')
1764  call hdf5_write(filehandle,crystallite_lp, 'Lp')
1765  call hdf5_write(filehandle,crystallite_li, 'Li')
1766  call hdf5_write(filehandle,crystallite_s, 'S')
1767 
1768  grouphandle = hdf5_addgroup(filehandle,'constituent')
1769  do i = 1,size(phase_plasticity)
1770  write(datasetname,'(i0,a)') i,'_omega_plastic'
1771  call hdf5_write(grouphandle,plasticstate(i)%state,datasetname)
1772  enddo
1773  call hdf5_closegroup(grouphandle)
1774 
1775  grouphandle = hdf5_addgroup(filehandle,'materialpoint')
1776  do i = 1, material_nhomogenization
1777  write(datasetname,'(i0,a)') i,'_omega_homogenization'
1778  call hdf5_write(grouphandle,homogstate(i)%state,datasetname)
1779  enddo
1780  call hdf5_closegroup(grouphandle)
1781 
1782  call hdf5_closefile(filehandle)
1783 
1784 end subroutine crystallite_restartwrite
1785 
1786 
1787 !--------------------------------------------------------------------------------------------------
1789 ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
1790 !--------------------------------------------------------------------------------------------------
1791 subroutine crystallite_restartread
1792 
1793  integer :: i
1794  integer(HID_T) :: filehandle, grouphandle
1795  character(len=pStringLen) :: filename, datasetname
1796 
1797  write(6,'(/,a,i0,a)') ' reading restart information of increment from file'
1798 
1799  write(filename,'(a,i0,a)') trim(getsolverjobname())//'_',worldrank,'.hdf5'
1800  filehandle = hdf5_openfile(filename)
1801 
1802  call hdf5_read(filehandle,crystallite_f0, 'F')
1803  call hdf5_read(filehandle,crystallite_fp0,'Fp')
1804  call hdf5_read(filehandle,crystallite_fi0,'Fi')
1805  call hdf5_read(filehandle,crystallite_lp0,'Lp')
1806  call hdf5_read(filehandle,crystallite_li0,'Li')
1807  call hdf5_read(filehandle,crystallite_s0, 'S')
1808 
1809  grouphandle = hdf5_opengroup(filehandle,'constituent')
1810  do i = 1,size(phase_plasticity)
1811  write(datasetname,'(i0,a)') i,'_omega_plastic'
1812  call hdf5_read(grouphandle,plasticstate(i)%state0,datasetname)
1813  enddo
1814  call hdf5_closegroup(grouphandle)
1815 
1816  grouphandle = hdf5_opengroup(filehandle,'materialpoint')
1817  do i = 1, material_nhomogenization
1818  write(datasetname,'(i0,a)') i,'_omega_homogenization'
1819  call hdf5_read(grouphandle,homogstate(i)%state0,datasetname)
1820  enddo
1821  call hdf5_closegroup(grouphandle)
1822 
1823  call hdf5_closefile(filehandle)
1824 
1825 end subroutine crystallite_restartread
1826 
1827 
1828 !--------------------------------------------------------------------------------------------------
1830 ! ToDo: Any guessing for the current states possible?
1831 !--------------------------------------------------------------------------------------------------
1832 subroutine crystallite_forward
1833 
1834  integer :: i, j
1835 
1836  crystallite_f0 = crystallite_partionedf
1837  crystallite_fp0 = crystallite_fp
1838  crystallite_lp0 = crystallite_lp
1839  crystallite_fi0 = crystallite_fi
1840  crystallite_li0 = crystallite_li
1841  crystallite_s0 = crystallite_s
1842 
1843  do i = 1, size(plasticstate)
1844  plasticstate(i)%state0 = plasticstate(i)%state
1845  enddo
1846  do i = 1, size(sourcestate)
1847  do j = 1,phase_nsources(i)
1848  sourcestate(i)%p(j)%state0 = sourcestate(i)%p(j)%state
1849  enddo; enddo
1850  do i = 1, material_nhomogenization
1851  homogstate(i)%state0 = homogstate(i)%state
1852  thermalstate(i)%state0 = thermalstate(i)%state
1853  damagestate(i)%state0 = damagestate(i)%state
1854  enddo
1855 
1856 end subroutine crystallite_forward
1857 
1858 end module crystallite
1859 
material::sourcestate
type(tsourcestate), dimension(:), allocatable, public sourcestate
Definition: material.f90:139
crystallite::crystallite_init
subroutine, public crystallite_init
allocates and initialize per grain variables
Definition: crystallite.f90:125
fesolving::fesolving_execip
integer, dimension(2) fesolving_execip
for ping-pong scheme always range to max IP, otherwise one specific IP
Definition: FEsolving.f90:18
crystallite::crystallite_s0
real(preal), dimension(:,:,:,:,:), allocatable, public, protected crystallite_s0
2nd Piola-Kirchhoff stress vector at start of FE inc
Definition: crystallite.f90:43
crystallite::nonlocalconvergencecheck
subroutine nonlocalconvergencecheck
sets convergence flag for nonlocal calculations
Definition: crystallite.f90:1675
constitutive::constitutive_lpanditstangents
subroutine, public constitutive_lpanditstangents(Lp, dLp_dS, dLp_dFi, S, Fi, ipc, ip, el)
contains the constitutive equation for calculating the velocity gradient
Definition: constitutive.f90:462
material::material_phaseat
integer, dimension(:,:), allocatable, public, protected material_phaseat
phase ID of each element
Definition: material.f90:132
crystallite::integratestatefpi
subroutine integratestatefpi
integrate stress, state with adaptive 1st order explicit Euler method using Fixed Point Iteration to ...
Definition: crystallite.f90:979
crystallite::crystallite_dt
real(preal), dimension(:,:,:), allocatable, public crystallite_dt
requested time increment of each grain
Definition: crystallite.f90:35
numerics::numerics_integrator
integer, public, protected numerics_integrator
method used for state integration Default 1: fix-point iteration
Definition: numerics.f90:1470
crystallite::crystallite_subfi0
real(preal), dimension(:,:,:,:,:), allocatable crystallite_subfi0
intermediate def grad at start of crystallite inc
Definition: crystallite.f90:65
rotations
rotation storage and conversion
Definition: rotations.f90:53
constitutive::constitutive_collectdotstate
subroutine, public constitutive_collectdotstate(S, FArray, Fi, FpArray, subdt, ipc, ip, el)
contains the constitutive equation for calculating the rate of change of microstructure
Definition: constitutive.f90:717
crystallite::tnumerics
Definition: crystallite.f90:87
lattice::lattice_bcc_id
@, public lattice_bcc_id
Definition: lattice.f90:395
crystallite::crystallite_stress
logical function, dimension(discretization_nip, discretization_nelem), public crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
calculate stress (P)
Definition: crystallite.f90:284
crystallite::crystallite_restartwrite
subroutine, public crystallite_restartwrite
Write current restart information (Field and constitutive data) to file.
Definition: crystallite.f90:1757
prec::tol_math_check
real(preal), parameter tol_math_check
tolerance for internal math self-checks (rotation)
Definition: prec.f90:30
crystallite::crystallite_fe
real(preal), dimension(:,:,:,:,:), allocatable, public, protected crystallite_fe
current "elastic" def grad (end of converged time step)
Definition: crystallite.f90:43
crystallite::integratestateeuler
subroutine integratestateeuler
integrate state with 1st order explicit Euler method
Definition: crystallite.f90:1145
crystallite::crystallite_forward
subroutine, public crystallite_forward
Forward data after successful increment.
Definition: crystallite.f90:1839
material::material_phasememberat
integer, dimension(:,:,:), allocatable, public, protected material_phasememberat
position of the element within its phase instance
Definition: material.f90:134
io::io_error
subroutine, public io_error(error_ID, el, ip, g, instance, ext_msg)
write error statements to standard out and terminate the Marc/spectral run with exit #9xxx
Definition: IO.f90:305
crystallite::crystallite_sublp0
real(preal), dimension(:,:,:,:,:), allocatable crystallite_sublp0
plastic velocity grad at start of crystallite inc
Definition: crystallite.f90:65
rotations::rotation
Definition: rotations.f90:63
crystallite::crystallite_subf
real(preal), dimension(:,:,:,:,:), allocatable crystallite_subf
def grad to be reached at end of crystallite inc
Definition: crystallite.f90:65
crystallite::crystallite_todo
logical, dimension(:,:,:), allocatable crystallite_todo
flag to indicate need for further computation
Definition: crystallite.f90:76
config::config_phase
type(tpartitionedstringlist), dimension(:), allocatable, public, protected config_phase
Definition: config.f90:23
crystallite::crystallite_partionedfp0
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_partionedfp0
plastic def grad at start of homog inc
Definition: crystallite.f90:52
crystallite::crystallite_restartread
subroutine, public crystallite_restartread
Read data for restart.
Definition: crystallite.f90:1798
math::math_99to3333
pure real(preal) function, dimension(3, 3, 3, 3) math_99to3333(m99)
convert 9x9 matrix into 3x3x3x3 matrix
Definition: math.f90:758
crystallite::crystallite_p
real(preal), dimension(:,:,:,:,:), allocatable, public, protected crystallite_p
1st Piola-Kirchhoff stress per grain
Definition: crystallite.f90:43
math::math_3333to99
pure real(preal) function, dimension(9, 9) math_3333to99(m3333)
convert 3x3x3x3 matrix into 9x9 matrix
Definition: math.f90:741
math::math_identity2nd
pure real(preal) function, dimension(d, d) math_identity2nd(d)
second rank identity tensor of specified dimension
Definition: math.f90:240
hdf5_utilities
Definition: HDF5_utilities.f90:11
constitutive::constitutive_dependentstate
subroutine, public constitutive_dependentstate(F, Fp, ipc, ip, el)
calls microstructure function of the different constitutive models
Definition: constitutive.f90:425
damper
real(preal) pure function damper(current, previous, previous2)
calculate the damping for correction of state and dot state
Definition: crystallite.f90:1122
crystallite::crystallite_orientation
type(rotation), dimension(:,:,:), allocatable crystallite_orientation
current orientation
Definition: crystallite.f90:41
config::config_name_phase
character(len=pstringlen), dimension(:), allocatable, public, protected config_name_phase
name of each phase
Definition: config.f90:34
crystallite::statejump
logical function statejump(ipc, ip, el)
calculates a jump in the state according to the current state and the current stress returns true,...
Definition: crystallite.f90:1704
math::math_33to9
pure real(preal) function, dimension(9) math_33to9(m33)
convert 3x3 matrix into vector 9
Definition: math.f90:650
crystallite::crystallite_fi0
real(preal), dimension(:,:,:,:,:), allocatable, public, protected crystallite_fi0
intermediate def grad at start of FE inc
Definition: crystallite.f90:43
crystallite::crystallite_fi
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_fi
current intermediate def grad (end of converged time step)
Definition: crystallite.f90:52
crystallite::crystallite_li
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_li
current intermediate velocitiy grad (end of converged time step)
Definition: crystallite.f90:52
crystallite::crystallite_s
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_s
current 2nd Piola-Kirchhoff stress vector (end of converged time step)
Definition: crystallite.f90:52
material
Parses material config file, either solverJobName.materialConfig or material.config.
Definition: material.f90:11
crystallite::toutput
new requested output (per phase)
Definition: crystallite.f90:81
crystallite::crystallite_converged
logical, dimension(:,:,:), allocatable crystallite_converged
convergence flag
Definition: crystallite.f90:76
material::homogenization_maxngrains
integer, public, protected homogenization_maxngrains
max number of grains in any USED homogenization
Definition: material.f90:110
material::material_orientation0
type(rotation), dimension(:,:,:), allocatable, public, protected material_orientation0
initial orientation of each grain,IP,element
Definition: material.f90:149
crystallite::output_constituent
type(toutput), dimension(:), allocatable output_constituent
Definition: crystallite.f90:85
crystallite::crystallite_fp0
real(preal), dimension(:,:,:,:,:), allocatable, public, protected crystallite_fp0
plastic def grad at start of FE inc
Definition: crystallite.f90:43
constitutive::constitutive_lianditstangents
subroutine, public constitutive_lianditstangents(Li, dLi_dS, dLi_dFi, S, Fi, ipc, ip, el)
contains the constitutive equation for calculating the velocity gradient
Definition: constitutive.f90:532
lattice::lattice_iso_id
@, public lattice_iso_id
Definition: lattice.f90:395
config
Reads in the material configuration from file.
Definition: config.f90:13
crystallite::crystallite_stresstangent
subroutine, public crystallite_stresstangent
calculate tangent (dPdF)
Definition: crystallite.f90:437
config::config_deallocate
subroutine, public config_deallocate(what)
deallocates the linked lists that store the content of the configuration files
Definition: config.f90:290
material::homogenization_ngrains
integer, dimension(:), allocatable, public, protected homogenization_ngrains
number of grains in each homogenization
Definition: material.f90:113
crystallite::num
type(tnumerics) num
Definition: crystallite.f90:103
crystallite::crystallite_localplasticity
logical, dimension(:,:,:), allocatable crystallite_localplasticity
indicates this grain to have purely local constitutive law
Definition: crystallite.f90:76
crystallite::crystallite_lp0
real(preal), dimension(:,:,:,:,:), allocatable, public, protected crystallite_lp0
plastic velocitiy grad at start of FE inc
Definition: crystallite.f90:43
crystallite::crystallite_dpdf
real(preal), dimension(:,:,:,:,:,:,:), allocatable, public, protected crystallite_dpdf
current individual dPdF per grain (end of converged time step)
Definition: crystallite.f90:72
crystallite::crystallite_partionedlp0
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_partionedlp0
plastic velocity grad at start of homog inc
Definition: crystallite.f90:52
crystallite::crystallite_subdt
real(preal), dimension(:,:,:), allocatable crystallite_subdt
substepped time increment of each grain
Definition: crystallite.f90:37
crystallite
crystallite state integration functions and reporting of results
Definition: crystallite.f90:15
prec
setting precision for real and int type
Definition: prec.f90:13
math::math_det33
real(preal) pure function math_det33(m)
determinant of a 3x3 matrix
Definition: math.f90:623
crystallite::crystallite_li0
real(preal), dimension(:,:,:,:,:), allocatable, public, protected crystallite_li0
intermediate velocitiy grad at start of FE inc
Definition: crystallite.f90:43
discretization
spatial discretization
Definition: discretization.f90:9
crystallite::crystallite_results
subroutine, public crystallite_results
writes crystallite results to HDF5 output file
Definition: crystallite.f90:625
select_rotations
type(rotation) function, dimension(:), allocatable select_rotations(dataset, instance)
select rotations for output
Definition: crystallite.f90:725
discretization::discretization_nip
integer, public, protected discretization_nip
Definition: discretization.f90:17
lattice::lattice_hex_id
@, public lattice_hex_id
Definition: lattice.f90:395
results::results_writedataset
Definition: results.f90:25
math::math_mul3333xx3333
pure real(preal) function, dimension(3, 3, 3, 3) math_mul3333xx3333(A, B)
matrix multiplication 3333x3333 = 3333 (ijkl,klmn)
Definition: math.f90:385
io::io_warning
subroutine, public io_warning(warning_ID, el, ip, g, ext_msg)
writes warning statement to standard out
Definition: IO.f90:535
material::phase_localplasticity
logical, dimension(:), allocatable, public, protected phase_localplasticity
flags phases with local constitutive law
Definition: material.f90:152
io
input/output functions, partly depending on chosen solver
Definition: IO.f90:12
material::phase_plasticityinstance
integer, dimension(:), allocatable, public, protected phase_plasticityinstance
instance of particular plasticity of each phase
Definition: material.f90:113
prec::preal
integer, parameter preal
number with 15 significant digits, up to 1e+-307 (typically 64 bit)
Definition: prec.f90:20
numerics::usepingpong
logical, public, protected usepingpong
Definition: numerics.f90:1483
debug
Reading in and interpretating the debugging settings for the various modules.
Definition: debug.f90:12
crystallite::integratestate
procedure(), pointer integratestate
Definition: crystallite.f90:105
crystallite::crystallite_partioneds0
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_partioneds0
2nd Piola-Kirchhoff stress vector at start of homog inc
Definition: crystallite.f90:52
math::math_invert33
pure subroutine math_invert33(InvA, DetA, error, A)
Cramer inversion of 3x3 matrix (subroutine)
Definition: math.f90:454
crystallite::crystallite_partionedfi0
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_partionedfi0
intermediate def grad at start of homog inc
Definition: crystallite.f90:52
results::results_addgroup
integer(hid_t) function, public results_addgroup(groupName)
adds a new group to the results file
Definition: results.f90:154
lattice::lattice_fcc_id
@, public lattice_fcc_id
Definition: lattice.f90:395
lattice::lattice_bct_id
@, public lattice_bct_id
Definition: lattice.f90:395
crystallite::crystallite_subli0
real(preal), dimension(:,:,:,:,:), allocatable crystallite_subli0
intermediate velocity grad at start of crystallite inc
Definition: crystallite.f90:65
constitutive
elasticity, plasticity, internal microstructure state
Definition: constitutive.f90:10
constitutive::constitutive_initialfi
pure real(preal) function, dimension(3, 3), public constitutive_initialfi(ipc, ip, el)
collects initial intermediate deformation gradient
Definition: constitutive.f90:610
lattice
contains lattice structure definitions including Schmid matrices for slip, twin, trans,
Definition: lattice.f90:13
crystallite::crystallite_subfp0
real(preal), dimension(:,:,:,:,:), allocatable crystallite_subfp0
plastic def grad at start of crystallite inc
Definition: crystallite.f90:65
math::math_invert
subroutine math_invert(InvA, error, A)
invert quadratic matrix of arbitrary dimension
Definition: math.f90:521
crystallite::crystallite_partionedf0
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_partionedf0
def grad at start of homog inc
Definition: crystallite.f90:52
crystallite::crystallite_fp
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_fp
current plastic def grad (end of converged time step)
Definition: crystallite.f90:52
crystallite::crystallite_substep
real(preal), dimension(:,:,:), allocatable crystallite_substep
size of next integration step
Definition: crystallite.f90:37
damask_interface
Interfacing between the 1-based solvers and the material subroutines provided by DAMASK.
Definition: DAMASK_interface.f90:22
results
Definition: results.f90:11
select_tensors
real(preal) function, dimension(:,:,:), allocatable select_tensors(dataset, instance)
select tensors for output
Definition: crystallite.f90:698
config::config_numerics
type(tpartitionedstringlist), public, protected config_numerics
Definition: config.f90:30
material::material_homogenizationat
integer, dimension(:), allocatable, public, protected material_homogenizationat
homogenization ID of each element (copy of discretization_homogenizationAt)
Definition: material.f90:128
material::plasticstate
type(tplasticstate), dimension(:), allocatable, public plasticstate
Definition: material.f90:137
math
Mathematical library, including random number generation and tensor representations.
Definition: math.f90:12
crystallite::crystallite_subfrac
real(preal), dimension(:,:,:), allocatable crystallite_subfrac
already calculated fraction of increment
Definition: crystallite.f90:37
crystallite::crystallite_f0
real(preal), dimension(:,:,:,:,:), allocatable, public, protected crystallite_f0
def grad at start of FE inc
Definition: crystallite.f90:43
crystallite::converged
logical pure function converged(residuum, state, atol)
determines whether a point is converged
Definition: crystallite.f90:1686
math::math_rotationalpart
real(preal) function, dimension(3, 3) math_rotationalpart(m)
rotational part from polar decomposition of 3x3 tensor
Definition: math.f90:962
crystallite::integratestaterkck45
subroutine integratestaterkck45
integrate stress, state with 5th order Runge-Kutta Cash-Karp method with adaptive step size (use 5th ...
Definition: crystallite.f90:1494
crystallite::crystallite_partionedli0
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_partionedli0
intermediate velocity grad at start of homog inc
Definition: crystallite.f90:52
math::math_9to33
pure real(preal) function, dimension(3, 3) math_9to33(v9)
convert 9 vector into 3x3 matrix
Definition: math.f90:667
discretization::discretization_nelem
integer, public, protected discretization_nelem
Definition: discretization.f90:17
crystallite::crystallite_subf0
real(preal), dimension(:,:,:,:,:), allocatable crystallite_subf0
def grad at start of crystallite inc
Definition: crystallite.f90:65
crystallite::integratestress
logical function integratestress(ipc, ip, el, timeFraction)
calculation of stress (P) with time integration based on a residuum in Lp and intermediate accelerati...
Definition: crystallite.f90:755
crystallite::crystallite_requested
logical, dimension(:,:,:), allocatable, public crystallite_requested
used by upper level (homogenization) to request crystallite calculation
Definition: crystallite.f90:74
crystallite::integratestaterk4
subroutine integratestaterk4
integrate stress, state with 4th order explicit Runge Kutta method
Definition: crystallite.f90:1339
results::results_closegroup
subroutine, public results_closegroup(group_id)
close a group
Definition: results.f90:166
crystallite::crystallite_partionedf
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_partionedf
def grad to be reached at end of homog inc
Definition: crystallite.f90:52
crystallite::crystallite_orientations
subroutine, public crystallite_orientations
calculates orientations
Definition: crystallite.f90:573
crystallite::integratestateadaptiveeuler
subroutine integratestateadaptiveeuler
integrate stress, state with 1st order Euler method with adaptive step size
Definition: crystallite.f90:1219
fesolving::fesolving_execelem
integer, dimension(2) fesolving_execelem
for ping-pong scheme always whole range, otherwise one specific element
Definition: FEsolving.f90:18
material::phase_nsources
integer, dimension(:), allocatable, public, protected phase_nsources
number of source mechanisms active in each phase
Definition: material.f90:113
crystallite::crystallite_lp
real(preal), dimension(:,:,:,:,:), allocatable, public crystallite_lp
current plastic velocitiy grad (end of converged time step)
Definition: crystallite.f90:52
lattice::lattice_structure
integer(kind(lattice_undefined_id)), dimension(:), allocatable, public, protected lattice_structure
Definition: lattice.f90:408
math::math_inv33
pure real(preal) function, dimension(3, 3) math_inv33(A)
Cramer inversion of 3x3 matrix (function)
Definition: math.f90:435
numerics
Managing of parameters related to numerics.
Definition: numerics.f90:10
fesolving
global variables for flow control
Definition: FEsolving.f90:10
constitutive::constitutive_sanditstangents
subroutine, public constitutive_sanditstangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to the elastic/intermediat...
Definition: constitutive.f90:645
lattice::lattice_ort_id
@, public lattice_ort_id
Definition: lattice.f90:395
crystallite::crystallite_push33toref
real(preal) function, dimension(3, 3), public crystallite_push33toref(ipc, ip, el, tensor33)
Map 2nd order tensor to reference config.
Definition: crystallite.f90:605
math::math_i3
real(preal), dimension(3, 3), parameter math_i3
3x3 Identity
Definition: math.f90:32