DAMASK with grid solvers  Revision: v2.0.3-2204-gdb1f2151
The Düsseldorf Advanced Material Simulation Kit with Grid Solvers
constitutive_plastic_nonlocal.f90
Go to the documentation of this file.
1 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
2 # 1 "<built-in>"
3 # 1 "<command-line>"
4 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
5 !--------------------------------------------------------------------------------------------------
10 !--------------------------------------------------------------------------------------------------
11 submodule(constitutive) plastic_nonlocal
12  use geometry_plastic_nonlocal, only: &
13  nipneighbors => geometry_plastic_nonlocal_nipneighbors, &
14  ipneighborhood => geometry_plastic_nonlocal_ipneighborhood, &
15  ipvolume => geometry_plastic_nonlocal_ipvolume0, &
16  iparea => geometry_plastic_nonlocal_iparea0, &
17  ipareanormal => geometry_plastic_nonlocal_ipareanormal0
18 
19  real(pReal), parameter :: &
20  kB = 1.38e-23_preal
21 
22  ! storage order of dislocation types
23  integer, dimension(8), parameter :: &
24  sgl = [1,2,3,4,5,6,7,8]
25  integer, dimension(5), parameter :: &
26  edg = [1,2,5,6,9], & !< edge
27  scr = [3,4,7,8,10]
28  integer, dimension(4), parameter :: &
29  mob = [1,2,3,4], & !< mobile
30  imm = [5,6,7,8]
31  integer, dimension(2), parameter :: &
32  dip = [9,10], & !< dipole
33  imm_edg = imm(1:2), &
34  imm_scr = imm(3:4)
35  integer, parameter :: &
36  mob_edg_pos = 1, & !< mobile edge positive
37  mob_edg_neg = 2, &
38  mob_scr_pos = 3, &
39  mob_scr_neg = 4
40 
41  ! BEGIN DEPRECATED
42  integer, dimension(:,:,:), allocatable :: &
43  iRhoU, & !< state indices for unblocked density
44  iV, & !< state indices for dislcation velocities
45  iD
46  !END DEPRECATED
47 
48  real(pReal), dimension(:,:,:,:,:,:), allocatable :: &
49  compatibility
50 
51  type :: tinitialparameters
52  real(pReal) :: &
53  rhoSglScatter, & !< standard deviation of scatter in initial dislocation density
54  rhoSglRandom, &
55  rhoSglRandomBinning
56  real(pReal), dimension(:), allocatable :: &
57  rhoSglEdgePos0, & !< initial edge_pos dislocation density
58  rhoSglEdgeNeg0, & !< initial edge_neg dislocation density
59  rhoSglScrewPos0, & !< initial screw_pos dislocation density
60  rhoSglScrewNeg0, & !< initial screw_neg dislocation density
61  rhoDipEdge0, & !< initial edge dipole dislocation density
62  rhoDipScrew0
63  integer, dimension(:) ,allocatable :: &
64  N_sl
65  end type tinitialparameters
66 
67  type :: tparameters
68  real(pReal) :: &
69  atomicVolume, & !< atomic volume
70  Dsd0, & !< prefactor for self-diffusion coefficient
71  selfDiffusionEnergy, & !< activation enthalpy for diffusion
72  atol_rho, & !< absolute tolerance for dislocation density in state integration
73  significantRho, & !< density considered significant
74  significantN, & !< number of dislocations considered significant
75  doublekinkwidth, & !< width of a doubkle kink in multiples of the burgers vector length b
76  solidSolutionEnergy, & !< activation energy for solid solution in J
77  solidSolutionSize, & !< solid solution obstacle size in multiples of the burgers vector length
78  solidSolutionConcentration, & !< concentration of solid solution in atomic parts
79  p, & !< parameter for kinetic law (Kocks,Argon,Ashby)
80  q, & !< parameter for kinetic law (Kocks,Argon,Ashby)
81  viscosity, & !< viscosity for dislocation glide in Pa s
82  fattack, & !< attack frequency in Hz
83  surfaceTransmissivity, & !< transmissivity at free surface
84  grainboundaryTransmissivity, & !< transmissivity at grain boundary (identified by different texture)
85  CFLfactor, & !< safety factor for CFL flux condition
86  fEdgeMultiplication, & !< factor that determines how much edge dislocations contribute to multiplication (0...1)
87  linetensionEffect, &
88  edgeJogFactor, &
89  mu, &
90  nu
91  real(pReal), dimension(:), allocatable :: &
92  minDipoleHeight_edge, & !< minimum stable edge dipole height
93  minDipoleHeight_screw, & !< minimum stable screw dipole height
94  peierlsstress_edge, &
95  peierlsstress_screw, &
96  lambda0, & !< mean free path prefactor for each
97  burgers
98  real(pReal), dimension(:,:), allocatable :: &
99  slip_normal, &
100  slip_direction, &
101  slip_transverse, &
102  minDipoleHeight, & ! edge and screw
103  peierlsstress, & ! edge and screw
104  interactionSlipSlip ,& !< coefficients for slip-slip interaction
105  forestProjection_Edge, & !< matrix of forest projections of edge dislocations
106  forestProjection_Screw
107  real(pReal), dimension(:,:,:), allocatable :: &
108  Schmid, & !< Schmid contribution
109  nonSchmid_pos, &
110  nonSchmid_neg
111  integer :: &
112  sum_N_sl
113  integer, dimension(:), allocatable :: &
114  colinearSystem
115  character(len=pStringLen), dimension(:), allocatable :: &
116  output
117  logical :: &
118  shortRangeStressCorrection, & !< use of short range stress correction by excess density gradient term
119  nonSchmidActive = .false.
120  end type tparameters
121 
122  type :: tnonlocalmicrostructure
123  real(pReal), allocatable, dimension(:,:) :: &
124  tau_pass, &
125  tau_Back
126  end type tnonlocalmicrostructure
127 
128  type :: tnonlocalstate
129  real(pReal), pointer, dimension(:,:) :: &
130  rho, & ! < all dislocations
131  rhoSgl, &
132  rhoSglMobile, & ! iRhoU
133  rho_sgl_mob_edg_pos, &
134  rho_sgl_mob_edg_neg, &
135  rho_sgl_mob_scr_pos, &
136  rho_sgl_mob_scr_neg, &
137  rhoSglImmobile, &
138  rho_sgl_imm_edg_pos, &
139  rho_sgl_imm_edg_neg, &
140  rho_sgl_imm_scr_pos, &
141  rho_sgl_imm_scr_neg, &
142  rhoDip, &
143  rho_dip_edg, &
144  rho_dip_scr, &
145  rho_forest, &
146  gamma, &
147  v, &
148  v_edg_pos, &
149  v_edg_neg, &
150  v_scr_pos, &
151  v_scr_neg
152  end type tnonlocalstate
153 
154  type(tNonlocalState), allocatable, dimension(:) :: &
155  deltaState, &
156  dotState, &
157  state, &
158  state0
159 
160  type(tParameters), dimension(:), allocatable :: param
161 
162  type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure
163 
164 contains
165 
166 !--------------------------------------------------------------------------------------------------
169 !--------------------------------------------------------------------------------------------------
170 module subroutine plastic_nonlocal_init
171 
172  integer :: &
173  Ninstance, &
174  p, &
175  NipcMyPhase, &
176  sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
177  s1, s2, &
178  s, t, l
179  real(pReal), dimension(:), allocatable :: &
180  a
181  character(len=pStringLen) :: &
182  extmsg = ''
183  type(tInitialParameters) :: &
184  ini
185 
186  write(6,'(/,a)') ' <<<+- constitutive_'//plasticity_nonlocal_label//' init -+>>>'; flush(6)
187 
188  write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014'
189  write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2014.03.012'
190 
191  write(6,'(/,a)') ' Kords, Dissertation RWTH Aachen, 2014'
192  write(6,'(a)') ' http://publications.rwth-aachen.de/record/229993'
193 
194  ninstance = count(phase_plasticity == plasticity_nonlocal_id)
195  if (iand(debug_level(debug_constitutive),debug_levelbasic) /= 0) &
196  write(6,'(a16,1x,i5,/)') '# instances:',ninstance
197 
198  allocate(param(ninstance))
199  allocate(state(ninstance))
200  allocate(state0(ninstance))
201  allocate(dotstate(ninstance))
202  allocate(deltastate(ninstance))
203  allocate(microstructure(ninstance))
204 
205  do p=1, size(config_phase)
206  if (phase_plasticity(p) /= plasticity_nonlocal_id) cycle
207 
208  associate(prm => param(phase_plasticityinstance(p)), &
209  dot => dotstate(phase_plasticityinstance(p)), &
210  stt => state(phase_plasticityinstance(p)), &
211  st0 => state0(phase_plasticityinstance(p)), &
212  del => deltastate(phase_plasticityinstance(p)), &
213  dst => microstructure(phase_plasticityinstance(p)), &
214  config => config_phase(p))
215 
216  prm%output = config%getStrings('(output)',defaultval=emptystringarray)
217 
218  prm%atol_rho = config%getFloat('atol_rho',defaultval=1.0e4_preal)
219 
220  ! This data is read in already in lattice
221  prm%mu = lattice_mu(p)
222  prm%nu = lattice_nu(p)
223 
224  ini%N_sl = config%getInts('nslip',defaultval=emptyintarray)
225  prm%sum_N_sl = sum(abs(ini%N_sl))
226  slipactive: if (prm%sum_N_sl > 0) then
227  prm%Schmid = lattice_schmidmatrix_slip(ini%N_sl,config%getString('lattice_structure'),&
228  config%getFloat('c/a',defaultval=0.0_preal))
229 
230  if(trim(config%getString('lattice_structure')) == 'bcc') then
231  a = config%getFloats('nonschmid_coefficients',defaultval = emptyrealarray)
232  if(size(a) > 0) prm%nonSchmidActive = .true.
233  prm%nonSchmid_pos = lattice_nonschmidmatrix(ini%N_sl,a,+1)
234  prm%nonSchmid_neg = lattice_nonschmidmatrix(ini%N_sl,a,-1)
235  else
236  prm%nonSchmid_pos = prm%Schmid
237  prm%nonSchmid_neg = prm%Schmid
238  endif
239 
240  prm%interactionSlipSlip = lattice_interaction_slipbyslip(ini%N_sl, &
241  config%getFloats('interaction_slipslip'), &
242  config%getString('lattice_structure'))
243 
244  prm%forestProjection_edge = lattice_forestprojection_edge(ini%N_sl,config%getString('lattice_structure'),&
245  config%getFloat('c/a',defaultval=0.0_preal))
246  prm%forestProjection_screw = lattice_forestprojection_screw(ini%N_sl,config%getString('lattice_structure'),&
247  config%getFloat('c/a',defaultval=0.0_preal))
248 
249  prm%slip_direction = lattice_slip_direction(ini%N_sl,config%getString('lattice_structure'),&
250  config%getFloat('c/a',defaultval=0.0_preal))
251  prm%slip_transverse = lattice_slip_transverse(ini%N_sl,config%getString('lattice_structure'),&
252  config%getFloat('c/a',defaultval=0.0_preal))
253  prm%slip_normal = lattice_slip_normal(ini%N_sl,config%getString('lattice_structure'),&
254  config%getFloat('c/a',defaultval=0.0_preal))
255 
256  ! collinear systems (only for octahedral slip systems in fcc)
257  allocate(prm%colinearSystem(prm%sum_N_sl), source = -1)
258  do s1 = 1, prm%sum_N_sl
259  do s2 = 1, prm%sum_N_sl
260  if (all(deq0(math_cross(prm%slip_direction(1:3,s1),prm%slip_direction(1:3,s2)))) .and. &
261  any(dneq0(math_cross(prm%slip_normal (1:3,s1),prm%slip_normal (1:3,s2))))) &
262  prm%colinearSystem(s1) = s2
263  enddo
264  enddo
265 
266  ini%rhoSglEdgePos0 = config%getFloats('rhosgledgepos0', requiredsize=size(ini%N_sl))
267  ini%rhoSglEdgeNeg0 = config%getFloats('rhosgledgeneg0', requiredsize=size(ini%N_sl))
268  ini%rhoSglScrewPos0 = config%getFloats('rhosglscrewpos0', requiredsize=size(ini%N_sl))
269  ini%rhoSglScrewNeg0 = config%getFloats('rhosglscrewneg0', requiredsize=size(ini%N_sl))
270  ini%rhoDipEdge0 = config%getFloats('rhodipedge0', requiredsize=size(ini%N_sl))
271  ini%rhoDipScrew0 = config%getFloats('rhodipscrew0', requiredsize=size(ini%N_sl))
272 
273  prm%lambda0 = config%getFloats('lambda0', requiredsize=size(ini%N_sl))
274  prm%burgers = config%getFloats('burgers', requiredsize=size(ini%N_sl))
275 
276  prm%lambda0 = math_expand(prm%lambda0,ini%N_sl)
277  prm%burgers = math_expand(prm%burgers,ini%N_sl)
278 
279  prm%minDipoleHeight_edge = config%getFloats('minimumdipoleheightedge', requiredsize=size(ini%N_sl))
280  prm%minDipoleHeight_screw = config%getFloats('minimumdipoleheightscrew', requiredsize=size(ini%N_sl))
281  prm%minDipoleHeight_edge = math_expand(prm%minDipoleHeight_edge, ini%N_sl)
282  prm%minDipoleHeight_screw = math_expand(prm%minDipoleHeight_screw,ini%N_sl)
283  allocate(prm%minDipoleHeight(prm%sum_N_sl,2))
284  prm%minDipoleHeight(:,1) = prm%minDipoleHeight_edge
285  prm%minDipoleHeight(:,2) = prm%minDipoleHeight_screw
286 
287  prm%peierlsstress_edge = config%getFloats('peierlsstressedge', requiredsize=size(ini%N_sl))
288  prm%peierlsstress_screw = config%getFloats('peierlsstressscrew', requiredsize=size(ini%N_sl))
289  prm%peierlsstress_edge = math_expand(prm%peierlsstress_edge, ini%N_sl)
290  prm%peierlsstress_screw = math_expand(prm%peierlsstress_screw,ini%N_sl)
291  allocate(prm%peierlsstress(prm%sum_N_sl,2))
292  prm%peierlsstress(:,1) = prm%peierlsstress_edge
293  prm%peierlsstress(:,2) = prm%peierlsstress_screw
294 
295  prm%significantRho = config%getFloat('significantrho')
296  prm%significantN = config%getFloat('significantn', 0.0_preal)
297  prm%CFLfactor = config%getFloat('cflfactor',defaultval=2.0_preal)
298 
299  prm%atomicVolume = config%getFloat('atomicvolume')
300  prm%Dsd0 = config%getFloat('selfdiffusionprefactor') !,'dsd0'
301  prm%selfDiffusionEnergy = config%getFloat('selfdiffusionenergy') !,'qsd'
302  prm%linetensionEffect = config%getFloat('linetension')
303  prm%edgeJogFactor = config%getFloat('edgejog')!,'edgejogs'
304  prm%doublekinkwidth = config%getFloat('doublekinkwidth')
305  prm%solidSolutionEnergy = config%getFloat('solidsolutionenergy')
306  prm%solidSolutionSize = config%getFloat('solidsolutionsize')
307  prm%solidSolutionConcentration = config%getFloat('solidsolutionconcentration')
308 
309  prm%p = config%getFloat('p')
310  prm%q = config%getFloat('q')
311  prm%viscosity = config%getFloat('viscosity')
312  prm%fattack = config%getFloat('attackfrequency')
313 
314  ! ToDo: discuss logic
315  ini%rhoSglScatter = config%getFloat('rhosglscatter')
316  ini%rhoSglRandom = config%getFloat('rhosglrandom',0.0_preal)
317  if (config%keyExists('/rhosglrandom/')) &
318  ini%rhoSglRandomBinning = config%getFloat('rhosglrandombinning',0.0_preal) !ToDo: useful default?
319  ! if (rhoSglRandom(instance) < 0.0_pReal) &
320  ! if (rhoSglRandomBinning(instance) <= 0.0_pReal) &
321 
322  prm%surfaceTransmissivity = config%getFloat('surfacetransmissivity',defaultval=1.0_preal)
323  prm%grainboundaryTransmissivity = config%getFloat('grainboundarytransmissivity',defaultval=-1.0_preal)
324  prm%fEdgeMultiplication = config%getFloat('edgemultiplication')
325  prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/')
326 
327 !--------------------------------------------------------------------------------------------------
328 ! sanity checks
329  if (any(prm%burgers < 0.0_preal)) extmsg = trim(extmsg)//' burgers'
330  if (any(prm%lambda0 <= 0.0_preal)) extmsg = trim(extmsg)//' lambda0'
331 
332  if (any(ini%rhoSglEdgePos0 < 0.0_preal)) extmsg = trim(extmsg)//' rhoSglEdgePos0'
333  if (any(ini%rhoSglEdgeNeg0 < 0.0_preal)) extmsg = trim(extmsg)//' rhoSglEdgeNeg0'
334  if (any(ini%rhoSglScrewPos0 < 0.0_preal)) extmsg = trim(extmsg)//' rhoSglScrewPos0'
335  if (any(ini%rhoSglScrewNeg0 < 0.0_preal)) extmsg = trim(extmsg)//' rhoSglScrewNeg0'
336  if (any(ini%rhoDipEdge0 < 0.0_preal)) extmsg = trim(extmsg)//' rhoDipEdge0'
337  if (any(ini%rhoDipScrew0 < 0.0_preal)) extmsg = trim(extmsg)//' rhoDipScrew0'
338 
339  if (any(prm%peierlsstress < 0.0_preal)) extmsg = trim(extmsg)//' peierlsstress'
340  if (any(prm%minDipoleHeight < 0.0_preal)) extmsg = trim(extmsg)//' minDipoleHeight'
341 
342  if (prm%viscosity <= 0.0_preal) extmsg = trim(extmsg)//' viscosity'
343  if (prm%selfDiffusionEnergy <= 0.0_preal) extmsg = trim(extmsg)//' selfDiffusionEnergy'
344  if (prm%fattack <= 0.0_preal) extmsg = trim(extmsg)//' fattack'
345  if (prm%doublekinkwidth <= 0.0_preal) extmsg = trim(extmsg)//' doublekinkwidth'
346  if (prm%Dsd0 < 0.0_preal) extmsg = trim(extmsg)//' Dsd0'
347  if (prm%atomicVolume <= 0.0_preal) extmsg = trim(extmsg)//' atomicVolume' ! ToDo: in disloUCLA, the atomic volume is given as a factor
348 
349  if (prm%significantN < 0.0_preal) extmsg = trim(extmsg)//' significantN'
350  if (prm%significantrho < 0.0_preal) extmsg = trim(extmsg)//' significantrho'
351  if (prm%atol_rho < 0.0_preal) extmsg = trim(extmsg)//' atol_rho'
352  if (prm%CFLfactor < 0.0_preal) extmsg = trim(extmsg)//' CFLfactor'
353 
354  if (prm%p <= 0.0_preal .or. prm%p > 1.0_preal) extmsg = trim(extmsg)//' p'
355  if (prm%q < 1.0_preal .or. prm%q > 2.0_preal) extmsg = trim(extmsg)//' q'
356 
357  if (prm%linetensionEffect < 0.0_preal .or. prm%linetensionEffect > 1.0_preal) &
358  extmsg = trim(extmsg)//' linetensionEffect'
359  if (prm%edgeJogFactor < 0.0_preal .or. prm%edgeJogFactor > 1.0_preal) &
360  extmsg = trim(extmsg)//' edgeJogFactor'
361 
362  if (prm%solidSolutionEnergy <= 0.0_preal) extmsg = trim(extmsg)//' solidSolutionEnergy'
363  if (prm%solidSolutionSize <= 0.0_preal) extmsg = trim(extmsg)//' solidSolutionSize'
364  if (prm%solidSolutionConcentration <= 0.0_preal) extmsg = trim(extmsg)//' solidSolutionConcentration'
365 
366  if (prm%grainboundaryTransmissivity > 1.0_preal) extmsg = trim(extmsg)//' grainboundaryTransmissivity'
367  if (prm%surfaceTransmissivity < 0.0_preal .or. prm%surfaceTransmissivity > 1.0_preal) &
368  extmsg = trim(extmsg)//' surfaceTransmissivity'
369 
370  if (prm%fEdgeMultiplication < 0.0_preal .or. prm%fEdgeMultiplication > 1.0_preal) &
371  extmsg = trim(extmsg)//' fEdgeMultiplication'
372 
373  endif slipactive
374 
375 !--------------------------------------------------------------------------------------------------
376 ! allocate state arrays
377  nipcmyphase = count(material_phaseat==p) * discretization_nip
378  sizedotstate = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
379  'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
380  'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
381  'rhoSglScrewPosImmobile','rhoSglScrewNegImmobile', &
382  'rhoDipEdge ','rhoDipScrew ', &
383  'gamma ' ]) * prm%sum_N_sl
384  sizedependentstate = size([ 'rhoForest ']) * prm%sum_N_sl
385  sizestate = sizedotstate + sizedependentstate &
386  + size([ 'velocityEdgePos ','velocityEdgeNeg ', &
387  'velocityScrewPos ','velocityScrewNeg ', &
388  'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl
389  sizedeltastate = sizedotstate
390 
391  call material_allocateplasticstate(p,nipcmyphase,sizestate,sizedotstate,sizedeltastate)
392 
393  plasticstate(p)%nonlocal = .true.
394  plasticstate(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
395 
396  st0%rho => plasticstate(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
397  stt%rho => plasticstate(p)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
398  dot%rho => plasticstate(p)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
399  del%rho => plasticstate(p)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
400  plasticstate(p)%atol(1:10*prm%sum_N_sl) = prm%atol_rho
401 
402  stt%rhoSgl => plasticstate(p)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
403  dot%rhoSgl => plasticstate(p)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
404  del%rhoSgl => plasticstate(p)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
405 
406  stt%rhoSglMobile => plasticstate(p)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
407  dot%rhoSglMobile => plasticstate(p)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
408  del%rhoSglMobile => plasticstate(p)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
409 
410  stt%rho_sgl_mob_edg_pos => plasticstate(p)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
411  dot%rho_sgl_mob_edg_pos => plasticstate(p)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
412  del%rho_sgl_mob_edg_pos => plasticstate(p)%deltaState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
413 
414  stt%rho_sgl_mob_edg_neg => plasticstate(p)%state (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
415  dot%rho_sgl_mob_edg_neg => plasticstate(p)%dotState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
416  del%rho_sgl_mob_edg_neg => plasticstate(p)%deltaState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
417 
418  stt%rho_sgl_mob_scr_pos => plasticstate(p)%state (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
419  dot%rho_sgl_mob_scr_pos => plasticstate(p)%dotState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
420  del%rho_sgl_mob_scr_pos => plasticstate(p)%deltaState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
421 
422  stt%rho_sgl_mob_scr_neg => plasticstate(p)%state (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
423  dot%rho_sgl_mob_scr_neg => plasticstate(p)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
424  del%rho_sgl_mob_scr_neg => plasticstate(p)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
425 
426  stt%rhoSglImmobile => plasticstate(p)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
427  dot%rhoSglImmobile => plasticstate(p)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
428  del%rhoSglImmobile => plasticstate(p)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
429 
430  stt%rho_sgl_imm_edg_pos => plasticstate(p)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
431  dot%rho_sgl_imm_edg_pos => plasticstate(p)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
432  del%rho_sgl_imm_edg_pos => plasticstate(p)%deltaState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
433 
434  stt%rho_sgl_imm_edg_neg => plasticstate(p)%state (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
435  dot%rho_sgl_imm_edg_neg => plasticstate(p)%dotState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
436  del%rho_sgl_imm_edg_neg => plasticstate(p)%deltaState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
437 
438  stt%rho_sgl_imm_scr_pos => plasticstate(p)%state (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
439  dot%rho_sgl_imm_scr_pos => plasticstate(p)%dotState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
440  del%rho_sgl_imm_scr_pos => plasticstate(p)%deltaState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
441 
442  stt%rho_sgl_imm_scr_neg => plasticstate(p)%state (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
443  dot%rho_sgl_imm_scr_neg => plasticstate(p)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
444  del%rho_sgl_imm_scr_neg => plasticstate(p)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
445 
446  stt%rhoDip => plasticstate(p)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
447  dot%rhoDip => plasticstate(p)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
448  del%rhoDip => plasticstate(p)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
449 
450  stt%rho_dip_edg => plasticstate(p)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
451  dot%rho_dip_edg => plasticstate(p)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
452  del%rho_dip_edg => plasticstate(p)%deltaState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
453 
454  stt%rho_dip_scr => plasticstate(p)%state (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
455  dot%rho_dip_scr => plasticstate(p)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
456  del%rho_dip_scr => plasticstate(p)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
457 
458  stt%gamma => plasticstate(p)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:nipcmyphase)
459  dot%gamma => plasticstate(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:nipcmyphase)
460  del%gamma => plasticstate(p)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:nipcmyphase)
461  plasticstate(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = config%getFloat('atol_gamma', defaultval = 1.0e-2_preal)
462  if(any(plasticstate(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_preal)) &
463  extmsg = trim(extmsg)//' atol_gamma'
464  plasticstate(p)%slipRate => plasticstate(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:nipcmyphase)
465 
466  stt%rho_forest => plasticstate(p)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:nipcmyphase)
467  stt%v => plasticstate(p)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:nipcmyphase)
468  stt%v_edg_pos => plasticstate(p)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:nipcmyphase)
469  stt%v_edg_neg => plasticstate(p)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:nipcmyphase)
470  stt%v_scr_pos => plasticstate(p)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:nipcmyphase)
471  stt%v_scr_neg => plasticstate(p)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:nipcmyphase)
472 
473  allocate(dst%tau_pass(prm%sum_N_sl,nipcmyphase),source=0.0_preal)
474  allocate(dst%tau_back(prm%sum_N_sl,nipcmyphase),source=0.0_preal)
475  end associate
476 
477  if (nipcmyphase > 0) call stateinit(ini,p,nipcmyphase)
478  plasticstate(p)%state0 = plasticstate(p)%state
479 
480 !--------------------------------------------------------------------------------------------------
481 ! exit if any parameter is out of range
482  if (extmsg /= '') call io_error(211,ext_msg=trim(extmsg)//'('//plasticity_nonlocal_label//')')
483 
484  enddo
485 
486  allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nipneighbors,&
487  discretization_nip,discretization_nelem), source=0.0_preal)
488 
489 ! BEGIN DEPRECATED----------------------------------------------------------------------------------
490  allocate(irhou(maxval(param%sum_N_sl),4,ninstance), source=0)
491  allocate(iv(maxval(param%sum_N_sl),4,ninstance), source=0)
492  allocate(id(maxval(param%sum_N_sl),2,ninstance), source=0)
493 
494  initializeinstances: do p = 1, size(phase_plasticity)
495  nipcmyphase = count(material_phaseat==p) * discretization_nip
496  myphase2: if (phase_plasticity(p) == plasticity_nonlocal_id) then
497  l = 0
498  do t = 1,4
499  do s = 1,param(phase_plasticityinstance(p))%sum_N_sl
500  l = l + 1
501  irhou(s,t,phase_plasticityinstance(p)) = l
502  enddo
503  enddo
504  l = l + (4+2+1+1)*param(phase_plasticityinstance(p))%sum_N_sl ! immobile(4), dipole(2), shear, forest
505  do t = 1,4
506  do s = 1,param(phase_plasticityinstance(p))%sum_N_sl
507  l = l + 1
508  iv(s,t,phase_plasticityinstance(p)) = l
509  enddo
510  enddo
511  do t = 1,2
512  do s = 1,param(phase_plasticityinstance(p))%sum_N_sl
513  l = l + 1
514  id(s,t,phase_plasticityinstance(p)) = l
515  enddo
516  enddo
517  if (id(param(phase_plasticityinstance(p))%sum_N_sl,2,phase_plasticityinstance(p)) /= plasticstate(p)%sizeState) &
518  call io_error(0, ext_msg = 'state indices not properly set ('//plasticity_nonlocal_label//')')
519  endif myphase2
520  enddo initializeinstances
521 
522 end subroutine plastic_nonlocal_init
523 
524 
525 !--------------------------------------------------------------------------------------------------
527 !--------------------------------------------------------------------------------------------------
528 module subroutine plastic_nonlocal_dependentstate(f, fp, instance, of, ip, el)
529 
530  real(pReal), dimension(3,3), intent(in) :: &
531  F, &
532  Fp
533  integer, intent(in) :: &
534  instance, &
535  of, &
536  ip, &
537  el
538 
539  integer :: &
540  no, & !< neighbor offset
541  neighbor_el, & ! element number of neighboring material point
542  neighbor_ip, & ! integration point of neighboring material point
543  neighbor_instance, & ! instance of this plasticity of neighboring material point
544  c, & ! index of dilsocation character (edge, screw)
545  s, & ! slip system index
546  dir, &
547  n
548  real(pReal) :: &
549  FVsize, &
550  nRealNeighbors ! number of really existing neighbors
551  integer, dimension(2) :: &
552  neighbors
553  real(pReal), dimension(2) :: &
554  rhoExcessGradient, &
555  rhoExcessGradient_over_rho, &
556  rhoTotal
557  real(pReal), dimension(3) :: &
558  rhoExcessDifferences, &
559  normal_latticeConf
560  real(pReal), dimension(3,3) :: &
561  invFe, & !< inverse of elastic deformation gradient
562  invFp, & !< inverse of plastic deformation gradient
563  connections, &
564  invConnections
565  real(pReal), dimension(3,nIPneighbors) :: &
566  connection_latticeConf
567  real(pReal), dimension(2,param(instance)%sum_N_sl) :: &
568  rhoExcess
569  real(pReal), dimension(param(instance)%sum_N_sl) :: &
570  rho_edg_delta, &
571  rho_scr_delta
572  real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
573  rho, &
574  rho0, &
575  rho_neighbor0
576  real(pReal), dimension(param(instance)%sum_N_sl,param(instance)%sum_N_sl) :: &
577  myInteractionMatrix ! corrected slip interaction matrix
578  real(pReal), dimension(param(instance)%sum_N_sl,nIPneighbors) :: &
579  rho_edg_delta_neighbor, &
580  rho_scr_delta_neighbor
581  real(pReal), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: &
582  neighbor_rhoExcess, & ! excess density at neighboring material point
583  neighbor_rhoTotal ! total density at neighboring material point
584  real(pReal), dimension(3,param(instance)%sum_N_sl,2) :: &
585  m ! direction of dislocation motion
586 
587  associate(prm => param(instance),dst => microstructure(instance), stt => state(instance))
588 
589  rho = getrho(instance,of,ip,el)
590 
591  stt%rho_forest(:,of) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
592  + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2))
593 
594 
595  ! coefficients are corrected for the line tension effect
596  ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals)
597  if (any(lattice_structure(material_phaseat(1,el)) == [lattice_bcc_id,lattice_fcc_id])) then
598  myinteractionmatrix = prm%interactionSlipSlip &
599  * spread(( 1.0_preal - prm%linetensionEffect &
600  + prm%linetensionEffect &
601  * log(0.35_preal * prm%burgers * sqrt(max(stt%rho_forest(:,of),prm%significantRho))) &
602  / log(0.35_preal * prm%burgers * 1e6_preal))** 2.0_preal,2,prm%sum_N_sl)
603  else
604  myinteractionmatrix = prm%interactionSlipSlip
605  endif
606 
607  dst%tau_pass(:,of) = prm%mu * prm%burgers &
608  * sqrt(matmul(myinteractionmatrix,sum(abs(rho),2)))
609 
610 !*** calculate the dislocation stress of the neighboring excess dislocation densities
611 !*** zero for material points of local plasticity
612 
613  !#################################################################################################
614  ! ToDo: MD: this is most likely only correct for F_i = I
615  !#################################################################################################
616 
617  rho0 = getrho0(instance,of,ip,el)
618  if (.not. phase_localplasticity(material_phaseat(1,el)) .and. prm%shortRangeStressCorrection) then
619  invfp = math_inv33(fp)
620  invfe = matmul(fp,math_inv33(f))
621 
622  rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
623  rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
624 
625  rhoexcess(1,:) = rho_edg_delta
626  rhoexcess(2,:) = rho_scr_delta
627 
628  fvsize = ipvolume(ip,el) ** (1.0_preal/3.0_preal)
629 
630  !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
631 
632  nrealneighbors = 0.0_preal
633  neighbor_rhototal = 0.0_preal
634  do n = 1,nipneighbors
635  neighbor_el = ipneighborhood(1,n,ip,el)
636  neighbor_ip = ipneighborhood(2,n,ip,el)
637  no = material_phasememberat(1,neighbor_ip,neighbor_el)
638  if (neighbor_el > 0 .and. neighbor_ip > 0) then
639  neighbor_instance = phase_plasticityinstance(material_phaseat(1,neighbor_el))
640  if (neighbor_instance == instance) then
641 
642  nrealneighbors = nrealneighbors + 1.0_preal
643  rho_neighbor0 = getrho0(instance,no,neighbor_ip,neighbor_el)
644 
645  rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg)
646  rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg)
647 
648  neighbor_rhototal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2)
649  neighbor_rhototal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2)
650 
651  connection_latticeconf(1:3,n) = matmul(invfe, discretization_ipcoords(1:3,neighbor_el+neighbor_ip-1) &
652  - discretization_ipcoords(1:3,el+neighbor_ip-1))
653  normal_latticeconf = matmul(transpose(invfp), ipareanormal(1:3,n,ip,el))
654  if (math_inner(normal_latticeconf,connection_latticeconf(1:3,n)) < 0.0_preal) & ! neighboring connection points in opposite direction to face normal: must be periodic image
655  connection_latticeconf(1:3,n) = normal_latticeconf * ipvolume(ip,el)/iparea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
656  else
657  ! local neighbor or different lattice structure or different constitution instance -> use central values instead
658  connection_latticeconf(1:3,n) = 0.0_preal
659  rho_edg_delta_neighbor(:,n) = rho_edg_delta
660  rho_scr_delta_neighbor(:,n) = rho_scr_delta
661  endif
662  else
663  ! free surface -> use central values instead
664  connection_latticeconf(1:3,n) = 0.0_preal
665  rho_edg_delta_neighbor(:,n) = rho_edg_delta
666  rho_scr_delta_neighbor(:,n) = rho_scr_delta
667  endif
668  enddo
669 
670  neighbor_rhoexcess(1,:,:) = rho_edg_delta_neighbor
671  neighbor_rhoexcess(2,:,:) = rho_scr_delta_neighbor
672 
673  !* loop through the slip systems and calculate the dislocation gradient by
674  !* 1. interpolation of the excess density in the neighorhood
675  !* 2. interpolation of the dead dislocation density in the central volume
676  m(1:3,:,1) = prm%slip_direction
677  m(1:3,:,2) = -prm%slip_transverse
678 
679  do s = 1,prm%sum_N_sl
680 
681  ! gradient from interpolation of neighboring excess density ...
682  do c = 1,2
683  do dir = 1,3
684  neighbors(1) = 2 * dir - 1
685  neighbors(2) = 2 * dir
686  connections(dir,1:3) = connection_latticeconf(1:3,neighbors(1)) &
687  - connection_latticeconf(1:3,neighbors(2))
688  rhoexcessdifferences(dir) = neighbor_rhoexcess(c,s,neighbors(1)) &
689  - neighbor_rhoexcess(c,s,neighbors(2))
690  enddo
691  invconnections = math_inv33(connections)
692  if (all(deq0(invconnections))) call io_error(-1,ext_msg='back stress calculation: inversion error')
693 
694  rhoexcessgradient(c) = math_inner(m(1:3,s,c), matmul(invconnections,rhoexcessdifferences))
695  enddo
696 
697  ! ... plus gradient from deads ...
698  rhoexcessgradient(1) = rhoexcessgradient(1) + sum(rho(s,imm_edg)) / fvsize
699  rhoexcessgradient(2) = rhoexcessgradient(2) + sum(rho(s,imm_scr)) / fvsize
700 
701  ! ... normalized with the total density ...
702  rhototal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhototal(1,s,:))) / (1.0_preal + nrealneighbors)
703  rhototal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhototal(2,s,:))) / (1.0_preal + nrealneighbors)
704 
705  rhoexcessgradient_over_rho = 0.0_preal
706  where(rhototal > 0.0_preal) rhoexcessgradient_over_rho = rhoexcessgradient / rhototal
707 
708  ! ... gives the local stress correction when multiplied with a factor
709  dst%tau_back(s,of) = - prm%mu * prm%burgers(s) / (2.0_preal * pi) &
710  * ( rhoexcessgradient_over_rho(1) / (1.0_preal - prm%nu) &
711  + rhoexcessgradient_over_rho(2))
712  enddo
713  endif
714 
715 # 721 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
716 
717  end associate
718 
719 end subroutine plastic_nonlocal_dependentstate
720 
721 
722 !--------------------------------------------------------------------------------------------------
724 !--------------------------------------------------------------------------------------------------
725 module subroutine plastic_nonlocal_lpanditstangent(lp,dlp_dmp, &
726  mp,temperature,instance,of,ip,el)
727  real(pReal), dimension(3,3), intent(out) :: &
728  Lp
729  real(pReal), dimension(3,3,3,3), intent(out) :: &
730  dLp_dMp
731  integer, intent(in) :: &
732  instance, &
733  of, &
734  ip, & !< current integration point
735  el
736  real(pReal), intent(in) :: &
737  Temperature
738 
739  real(pReal), dimension(3,3), intent(in) :: &
740  Mp
741 
742  integer :: &
743  ns, & !< short notation for the total number of active slip systems
744  i, &
745  j, &
746  k, &
747  l, &
748  t, & !< dislocation type
749  s
750  real(pReal), dimension(param(instance)%sum_N_sl,8) :: &
751  rhoSgl
752  real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
753  rho
754  real(pReal), dimension(param(instance)%sum_N_sl,4) :: &
755  v, & !< velocity
756  tauNS, & !< resolved shear stress including non Schmid and backstress terms
757  dv_dtau, & !< velocity derivative with respect to the shear stress
758  dv_dtauNS
759  real(pReal), dimension(param(instance)%sum_N_sl) :: &
760  tau, & !< resolved shear stress including backstress terms
761  gdotTotal
762 
763  associate(prm => param(instance),dst=>microstructure(instance),stt=>state(instance))
764  ns = prm%sum_N_sl
765 
766  !*** shortcut to state variables
767  rho = getrho(instance,of,ip,el)
768  rhosgl = rho(:,sgl)
769 
770  do s = 1,ns
771  tau(s) = math_tensordot(mp, prm%Schmid(1:3,1:3,s))
772  tauns(s,1) = tau(s)
773  tauns(s,2) = tau(s)
774  if (tau(s) > 0.0_preal) then
775  tauns(s,3) = math_tensordot(mp, +prm%nonSchmid_pos(1:3,1:3,s))
776  tauns(s,4) = math_tensordot(mp, -prm%nonSchmid_neg(1:3,1:3,s))
777  else
778  tauns(s,3) = math_tensordot(mp, +prm%nonSchmid_neg(1:3,1:3,s))
779  tauns(s,4) = math_tensordot(mp, -prm%nonSchmid_pos(1:3,1:3,s))
780  endif
781  enddo
782  tauns = tauns + spread(dst%tau_back(:,of),2,4)
783  tau = tau + dst%tau_back(:,of)
784 
785  ! edges
786  call kinetics(v(:,1), dv_dtau(:,1), dv_dtauns(:,1), &
787  tau, tauns(:,1), dst%tau_pass(:,of),1,temperature, instance)
788  v(:,2) = v(:,1)
789  dv_dtau(:,2) = dv_dtau(:,1)
790  dv_dtauns(:,2) = dv_dtauns(:,1)
791 
792  !screws
793  if (prm%nonSchmidActive) then
794  v(:,3:4) = spread(v(:,1),2,2)
795  dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2)
796  dv_dtauns(:,3:4) = spread(dv_dtauns(:,1),2,2)
797  else
798  do t = 3,4
799  call kinetics(v(:,t), dv_dtau(:,t), dv_dtauns(:,t), &
800  tau, tauns(:,t), dst%tau_pass(:,of),2,temperature, instance)
801  enddo
802  endif
803 
804  stt%v(:,of) = pack(v,.true.)
805 
806  !*** Bauschinger effect
807  forall (s = 1:ns, t = 5:8, rhosgl(s,t) * v(s,t-4) < 0.0_preal) &
808  rhosgl(s,t-4) = rhosgl(s,t-4) + abs(rhosgl(s,t))
809 
810  gdottotal = sum(rhosgl(:,1:4) * v, 2) * prm%burgers
811 
812  lp = 0.0_preal
813  dlp_dmp = 0.0_preal
814  do s = 1,ns
815  lp = lp + gdottotal(s) * prm%Schmid(1:3,1:3,s)
816  forall (i=1:3,j=1:3,k=1:3,l=1:3) &
817  dlp_dmp(i,j,k,l) = dlp_dmp(i,j,k,l) &
818  + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) &
819  * sum(rhosgl(s,1:4) * dv_dtau(s,1:4)) * prm%burgers(s) &
820  + prm%Schmid(i,j,s) &
821  * ( prm%nonSchmid_pos(k,l,s) * rhosgl(s,3) * dv_dtauns(s,3) &
822  - prm%nonSchmid_neg(k,l,s) * rhosgl(s,4) * dv_dtauns(s,4)) * prm%burgers(s)
823  enddo
824 
825  end associate
826 
827 end subroutine plastic_nonlocal_lpanditstangent
828 
829 
830 !--------------------------------------------------------------------------------------------------
832 !--------------------------------------------------------------------------------------------------
833 module subroutine plastic_nonlocal_deltastate(mp,instance,of,ip,el)
834 
835  real(pReal), dimension(3,3), intent(in) :: &
836  Mp
837  integer, intent(in) :: &
838  instance, & ! current instance of this plasticity
839  of, & !< offset
840  ip, &
841  el
842 
843  integer :: &
844  ph, & !< phase
845  ns, & ! short notation for the total number of active slip systems
846  c, & ! character of dislocation
847  t, & ! type of dislocation
848  s ! index of my current slip system
849  real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
850  deltaRhoRemobilization, & ! density increment by remobilization
851  deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change)
852  real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
853  rho ! current dislocation densities
854  real(pReal), dimension(param(instance)%sum_N_sl,4) :: &
855  v ! dislocation glide velocity
856  real(pReal), dimension(param(instance)%sum_N_sl) :: &
857  tau ! current resolved shear stress
858  real(pReal), dimension(param(instance)%sum_N_sl,2) :: &
859  rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
860  dUpper, & ! current maximum stable dipole distance for edges and screws
861  dUpperOld, & ! old maximum stable dipole distance for edges and screws
862  deltaDUpper ! change in maximum stable dipole distance for edges and screws
863 
864  ph = material_phaseat(1,el)
865 
866  associate(prm => param(instance),dst => microstructure(instance),del => deltastate(instance))
867  ns = prm%sum_N_sl
868 
869  !*** shortcut to state variables
870  forall (s = 1:ns, t = 1:4) v(s,t) = plasticstate(ph)%state(iv(s,t,instance),of)
871  forall (s = 1:ns, c = 1:2) dupperold(s,c) = plasticstate(ph)%state(id(s,c,instance),of)
872 
873  rho = getrho(instance,of,ip,el)
874  rhodip = rho(:,dip)
875 
876  !****************************************************************************
877  !*** dislocation remobilization (bauschinger effect)
878  where(rho(:,imm) * v < 0.0_preal)
879  deltarhoremobilization(:,mob) = abs(rho(:,imm))
880  deltarhoremobilization(:,imm) = - rho(:,imm)
881  rho(:,mob) = rho(:,mob) + abs(rho(:,imm))
882  rho(:,imm) = 0.0_preal
883  elsewhere
884  deltarhoremobilization(:,mob) = 0.0_preal
885  deltarhoremobilization(:,imm) = 0.0_preal
886  endwhere
887  deltarhoremobilization(:,dip) = 0.0_preal
888 
889  !****************************************************************************
890  !*** calculate dipole formation and dissociation by stress change
891 
892  !*** calculate limits for stable dipole height
893  do s = 1,prm%sum_N_sl
894  tau(s) = math_tensordot(mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,of)
895  if (abs(tau(s)) < 1.0e-15_preal) tau(s) = 1.0e-15_preal
896  enddo
897 
898  dupper(:,1) = prm%mu * prm%burgers/(8.0_preal * pi * (1.0_preal - prm%nu) * abs(tau))
899  dupper(:,2) = prm%mu * prm%burgers/(4.0_preal * pi * abs(tau))
900 
901  where(dneq0(sqrt(sum(abs(rho(:,edg)),2)))) &
902  dupper(:,1) = min(1.0_preal/sqrt(sum(abs(rho(:,edg)),2)),dupper(:,1))
903  where(dneq0(sqrt(sum(abs(rho(:,scr)),2)))) &
904  dupper(:,2) = min(1.0_preal/sqrt(sum(abs(rho(:,scr)),2)),dupper(:,2))
905 
906  dupper = max(dupper,prm%minDipoleHeight)
907  deltadupper = dupper - dupperold
908 
909 
910  !*** dissociation by stress increase
911  deltarhodipole2singlestress = 0.0_preal
912  forall (c=1:2, s=1:ns, deltadupper(s,c) < 0.0_preal .and. &
913  dneq0(dupperold(s,c) - prm%minDipoleHeight(s,c))) &
914  deltarhodipole2singlestress(s,8+c) = rhodip(s,c) * deltadupper(s,c) &
915  / (dupperold(s,c) - prm%minDipoleHeight(s,c))
916 
917  forall (t=1:4) deltarhodipole2singlestress(:,t) = -0.5_preal * deltarhodipole2singlestress(:,(t-1)/2+9)
918  forall (s = 1:ns, c = 1:2) plasticstate(ph)%state(id(s,c,instance),of) = dupper(s,c)
919 
920  plasticstate(ph)%deltaState(:,of) = 0.0_preal
921  del%rho(:,of) = reshape(deltarhoremobilization + deltarhodipole2singlestress, [10*ns])
922 
923 # 936 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
924 
925  end associate
926 
927 end subroutine plastic_nonlocal_deltastate
928 
929 
930 !---------------------------------------------------------------------------------------------------
932 !---------------------------------------------------------------------------------------------------
933 module subroutine plastic_nonlocal_dotstate(mp, f, fp, temperature,timestep, &
934  instance,of,ip,el)
935 
936  real(pReal), dimension(3,3), intent(in) :: &
937  Mp
938  real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: &
939  F, & !< elastic deformation gradient
940  Fp
941  real(pReal), intent(in) :: &
942  Temperature, & !< temperature
943  timestep
944  integer, intent(in) :: &
945  instance, &
946  of, &
947  ip, & !< current integration point
948  el
949 
950  integer :: &
951  ph, &
952  neighbor_instance, & !< instance of my neighbor's plasticity
953  ns, & !< short notation for the total number of active slip systems
954  c, & !< character of dislocation
955  n, & !< index of my current neighbor
956  neighbor_el, & !< element number of my neighbor
957  neighbor_ip, & !< integration point of my neighbor
958  neighbor_n, & !< neighbor index pointing to me when looking from my neighbor
959  opposite_neighbor, & !< index of my opposite neighbor
960  opposite_ip, & !< ip of my opposite neighbor
961  opposite_el, & !< element index of my opposite neighbor
962  opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor
963  t, & !< type of dislocation
964  no,& !< neighbor offset shortcut
965  np,& !< neighbor phase shortcut
966  topp, & !< type of dislocation with opposite sign to t
967  s !< index of my current slip system
968  real(pReal), dimension(param(instance)%sum_N_sl,10) :: &
969  rho, &
970  rho0, & !< dislocation density at beginning of time step
971  rhoDot, & !< density evolution
972  rhoDotMultiplication, & !< density evolution by multiplication
973  rhoDotFlux, & !< density evolution by flux
974  rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
975  rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
976  rhoDotThermalAnnihilation !< density evolution by thermal annihilation
977  real(pReal), dimension(param(instance)%sum_N_sl,8) :: &
978  rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
979  neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
980  my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
981  real(pReal), dimension(param(instance)%sum_N_sl,4) :: &
982  v, & !< current dislocation glide velocity
983  v0, &
984  neighbor_v0, & !< dislocation glide velocity of enighboring ip
985  gdot !< shear rates
986  real(pReal), dimension(param(instance)%sum_N_sl) :: &
987  tau, & !< current resolved shear stress
988  vClimb !< climb velocity of edge dipoles
989  real(pReal), dimension(param(instance)%sum_N_sl,2) :: &
990  rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
991  dLower, & !< minimum stable dipole distance for edges and screws
992  dUpper !< current maximum stable dipole distance for edges and screws
993  real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: &
994  m !< direction of dislocation motion
995  real(pReal), dimension(3,3) :: &
996  my_F, & !< my total deformation gradient
997  neighbor_F, & !< total deformation gradient of my neighbor
998  my_Fe, & !< my elastic deformation gradient
999  neighbor_Fe, & !< elastic deformation gradient of my neighbor
1000  Favg !< average total deformation gradient of me and my neighbor
1001  real(pReal), dimension(3) :: &
1002  normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration
1003  normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration
1004  normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration
1005  normal_me2neighbor_defConf
1006  real(pReal) :: &
1007  area, & !< area of the current interface
1008  transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point
1009  lineLength, & !< dislocation line length leaving the current interface
1010  selfDiffusion
1011 
1012  ph = material_phaseat(1,el)
1013  if (timestep <= 0.0_preal) then
1014  plasticstate(ph)%dotState = 0.0_preal
1015  return
1016  endif
1017 
1018  associate(prm => param(instance), &
1019  dst => microstructure(instance), &
1020  dot => dotstate(instance), &
1021  stt => state(instance))
1022  ns = prm%sum_N_sl
1023 
1024  tau = 0.0_preal
1025  gdot = 0.0_preal
1026 
1027  rho = getrho(instance,of,ip,el)
1028  rhosgl = rho(:,sgl)
1029  rhodip = rho(:,dip)
1030  rho0 = getrho0(instance,of,ip,el)
1031  my_rhosgl0 = rho0(:,sgl)
1032 
1033  forall (s = 1:ns, t = 1:4) v(s,t) = plasticstate(ph)%state(iv(s,t,instance),of)
1034  gdot = rhosgl(:,1:4) * v * spread(prm%burgers,2,4)
1035 
1036 # 1056 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
1037 
1038  !****************************************************************************
1039  !*** limits for stable dipole height
1040  do s = 1,ns
1041  tau(s) = math_tensordot(mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,of)
1042  if (abs(tau(s)) < 1.0e-15_preal) tau(s) = 1.0e-15_preal
1043  enddo
1044 
1045  dlower = prm%minDipoleHeight
1046  dupper(:,1) = prm%mu * prm%burgers/(8.0_preal * pi * (1.0_preal - prm%nu) * abs(tau))
1047  dupper(:,2) = prm%mu * prm%burgers/(4.0_preal * pi * abs(tau))
1048 
1049  where(dneq0(sqrt(sum(abs(rho(:,edg)),2)))) &
1050  dupper(:,1) = min(1.0_preal/sqrt(sum(abs(rho(:,edg)),2)),dupper(:,1))
1051  where(dneq0(sqrt(sum(abs(rho(:,scr)),2)))) &
1052  dupper(:,2) = min(1.0_preal/sqrt(sum(abs(rho(:,scr)),2)),dupper(:,2))
1053 
1054  dupper = max(dupper,dlower)
1055 
1056  !****************************************************************************
1057  !*** dislocation multiplication
1058  rhodotmultiplication = 0.0_preal
1059  isbcc: if (lattice_structure(ph) == lattice_bcc_id) then
1060  forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_preal)
1061  rhodotmultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
1062  * sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path
1063  ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation
1064  rhodotmultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
1065  * sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path
1066  ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation
1067  endforall
1068 
1069  else isbcc
1070  rhodotmultiplication(:,1:4) = spread( &
1071  (sum(abs(gdot(:,1:2)),2) * prm%fEdgeMultiplication + sum(abs(gdot(:,3:4)),2)) &
1072  * sqrt(stt%rho_forest(:,of)) / prm%lambda0 / prm%burgers, 2, 4)
1073  endif isbcc
1074 
1075  forall (s = 1:ns, t = 1:4) v0(s,t) = plasticstate(ph)%state0(iv(s,t,instance),of)
1076 
1077  !****************************************************************************
1078  !*** calculate dislocation fluxes (only for nonlocal plasticity)
1079  rhodotflux = 0.0_preal
1080  if (.not. phase_localplasticity(material_phaseat(1,el))) then
1081 
1082  !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
1083  if (any( abs(gdot) > 0.0_preal & ! any active slip system ...
1084  .and. prm%CFLfactor * abs(v0) * timestep &
1085  > ipvolume(ip,el) / maxval(iparea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
1086 # 1116 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
1087  plasticstate(ph)%dotState = ieee_value(1.0_preal,ieee_quiet_nan) ! -> return NaN and, hence, enforce cutback
1088  return
1089  endif
1090 
1091 
1092  !*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!!
1093  !*** opposite sign to our t vector in the (s,t,n) triplet !!!
1094 
1095  m(1:3,:,1) = prm%slip_direction
1096  m(1:3,:,2) = -prm%slip_direction
1097  m(1:3,:,3) = -prm%slip_transverse
1098  m(1:3,:,4) = prm%slip_transverse
1099 
1100  my_f = f(1:3,1:3,1,ip,el)
1101  my_fe = matmul(my_f, math_inv33(fp(1:3,1:3,1,ip,el)))
1102 
1103  neighbors: do n = 1,nipneighbors
1104 
1105  neighbor_el = ipneighborhood(1,n,ip,el)
1106  neighbor_ip = ipneighborhood(2,n,ip,el)
1107  neighbor_n = ipneighborhood(3,n,ip,el)
1108  np = material_phaseat(1,neighbor_el)
1109  no = material_phasememberat(1,neighbor_ip,neighbor_el)
1110 
1111  opposite_neighbor = n + mod(n,2) - mod(n+1,2)
1112  opposite_el = ipneighborhood(1,opposite_neighbor,ip,el)
1113  opposite_ip = ipneighborhood(2,opposite_neighbor,ip,el)
1114  opposite_n = ipneighborhood(3,opposite_neighbor,ip,el)
1115 
1116  if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
1117  neighbor_instance = phase_plasticityinstance(material_phaseat(1,neighbor_el))
1118  neighbor_f = f(1:3,1:3,1,neighbor_ip,neighbor_el)
1119  neighbor_fe = matmul(neighbor_f, math_inv33(fp(1:3,1:3,1,neighbor_ip,neighbor_el)))
1120  favg = 0.5_preal * (my_f + neighbor_f)
1121  else ! if no neighbor, take my value as average
1122  favg = my_f
1123  endif
1124 
1125  neighbor_v0 = 0.0_preal ! needed for check of sign change in flux density below
1126 
1127  !* FLUX FROM MY NEIGHBOR TO ME
1128  !* This is only considered, if I have a neighbor of nonlocal plasticity
1129  !* (also nonlocal constitutive law with local properties) that is at least a little bit
1130  !* compatible.
1131  !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of
1132  !* my neighbor's interface.
1133  !* The entering flux from my neighbor will be distributed on my slip systems according to the
1134  !* compatibility
1135  if (neighbor_n > 0) then
1136  if (phase_plasticity(material_phaseat(1,neighbor_el)) == plasticity_nonlocal_id .and. &
1137  any(compatibility(:,:,:,n,ip,el) > 0.0_preal)) then
1138 
1139  forall (s = 1:ns, t = 1:4)
1140  neighbor_v0(s,t) = plasticstate(np)%state0(iv(s,t,neighbor_instance),no)
1141  neighbor_rhosgl0(s,t) = max(plasticstate(np)%state0(irhou(s,t,neighbor_instance),no),0.0_preal)
1142  endforall
1143 
1144  where (neighbor_rhosgl0 * ipvolume(neighbor_ip,neighbor_el) ** 0.667_preal < prm%significantN &
1145  .or. neighbor_rhosgl0 < prm%significantRho) &
1146  neighbor_rhosgl0 = 0.0_preal
1147  normal_neighbor2me_defconf = math_det33(favg) * matmul(math_inv33(transpose(favg)), &
1148  ipareanormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => me)
1149  normal_neighbor2me = matmul(transpose(neighbor_fe), normal_neighbor2me_defconf) &
1150  / math_det33(neighbor_fe) ! interface normal in the lattice configuration of my neighbor
1151  area = iparea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me)
1152  normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length
1153  do s = 1,ns
1154  do t = 1,4
1155  c = (t + 1) / 2
1156  topp = t + mod(t,2) - mod(t+1,2)
1157  if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_preal & ! flux from my neighbor to me == entering flux for me
1158  .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_preal ) then ! ... only if no sign change in flux density
1159  linelength = neighbor_rhosgl0(s,t) * neighbor_v0(s,t) &
1160  * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
1161  where (compatibility(c,:,s,n,ip,el) > 0.0_preal) &
1162  rhodotflux(:,t) = rhodotflux(1:ns,t) &
1163  + linelength/ipvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_preal ! transferring to equally signed mobile dislocation type
1164  where (compatibility(c,:,s,n,ip,el) < 0.0_preal) &
1165  rhodotflux(:,topp) = rhodotflux(:,topp) &
1166  + linelength/ipvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_preal ! transferring to opposite signed mobile dislocation type
1167 
1168  endif
1169  enddo
1170  enddo
1171  endif; endif
1172 
1173 
1174  !* FLUX FROM ME TO MY NEIGHBOR
1175  !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties).
1176  !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me.
1177  !* So the net flux in the direction of my neighbor is equal to zero:
1178  !* leaving flux to neighbor == entering flux from opposite neighbor
1179  !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
1180  !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
1181  if (opposite_n > 0) then
1182  if (phase_plasticity(material_phaseat(1,opposite_el)) == plasticity_nonlocal_id) then
1183 
1184  normal_me2neighbor_defconf = math_det33(favg) &
1185  * matmul(math_inv33(transpose(favg)),ipareanormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing me => neighbor)
1186  normal_me2neighbor = matmul(transpose(my_fe), normal_me2neighbor_defconf) &
1187  / math_det33(my_fe) ! interface normal in my lattice configuration
1188  area = iparea(n,ip,el) * norm2(normal_me2neighbor)
1189  normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length
1190  do s = 1,ns
1191  do t = 1,4
1192  c = (t + 1) / 2
1193  if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_preal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
1194  if (v0(s,t) * neighbor_v0(s,t) >= 0.0_preal) then ! no sign change in flux density
1195  transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_preal) ! overall transmissivity from this slip system to my neighbor
1196  else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
1197  transmissivity = 0.0_preal
1198  endif
1199  linelength = my_rhosgl0(s,t) * v0(s,t) &
1200  * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
1201  rhodotflux(s,t) = rhodotflux(s,t) - linelength / ipvolume(ip,el) ! subtract dislocation flux from current type
1202  rhodotflux(s,t+4) = rhodotflux(s,t+4) &
1203  + linelength / ipvolume(ip,el) * (1.0_preal - transmissivity) &
1204  * sign(1.0_preal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
1205  endif
1206  enddo
1207  enddo
1208  endif; endif
1209 
1210  enddo neighbors
1211  endif
1212 
1213 
1214 
1215  !****************************************************************************
1216  !*** calculate dipole formation and annihilation
1217 
1218  !*** formation by glide
1219  do c = 1,2
1220  rhodotsingle2dipoleglide(:,2*c-1) = -2.0_preal * dupper(:,c) / prm%burgers &
1221  * ( rhosgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile
1222  + rhosgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile
1223  + abs(rhosgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile
1224 
1225  rhodotsingle2dipoleglide(:,2*c) = -2.0_preal * dupper(:,c) / prm%burgers &
1226  * ( rhosgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile
1227  + rhosgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile
1228  + abs(rhosgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile
1229 
1230  rhodotsingle2dipoleglide(:,2*c+3) = -2.0_preal * dupper(:,c) / prm%burgers &
1231  * rhosgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile
1232 
1233  rhodotsingle2dipoleglide(:,2*c+4) = -2.0_preal * dupper(:,c) / prm%burgers &
1234  * rhosgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile
1235 
1236  rhodotsingle2dipoleglide(:,c+8) = abs(rhodotsingle2dipoleglide(:,2*c+3)) &
1237  + abs(rhodotsingle2dipoleglide(:,2*c+4)) &
1238  - rhodotsingle2dipoleglide(:,2*c-1) &
1239  - rhodotsingle2dipoleglide(:,2*c)
1240  enddo
1241 
1242 
1243  !*** athermal annihilation
1244  rhodotathermalannihilation = 0.0_preal
1245  forall (c=1:2) &
1246  rhodotathermalannihilation(:,c+8) = -2.0_preal * dlower(:,c) / prm%burgers &
1247  * ( 2.0_preal * (rhosgl(:,2*c-1) * abs(gdot(:,2*c)) + rhosgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single
1248  + 2.0_preal * (abs(rhosgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhosgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
1249  + rhodip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent
1250 
1251  ! annihilated screw dipoles leave edge jogs behind on the colinear system
1252  if (lattice_structure(ph) == lattice_fcc_id) &
1253  forall (s = 1:ns, prm%colinearSystem(s) > 0) &
1254  rhodotathermalannihilation(prm%colinearSystem(s),1:2) = - rhodotathermalannihilation(s,10) &
1255  * 0.25_preal * sqrt(stt%rho_forest(s,of)) * (dupper(s,2) + dlower(s,2)) * prm%edgeJogFactor
1256 
1257 
1258  !*** thermally activated annihilation of edge dipoles by climb
1259  rhodotthermalannihilation = 0.0_preal
1260  selfdiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (kb * temperature))
1261  vclimb = prm%atomicVolume * selfdiffusion * prm%mu &
1262  / ( kb * temperature * pi * (1.0_preal-prm%nu) * (dupper(:,1) + dlower(:,1)))
1263  forall (s = 1:ns, dupper(s,1) > dlower(s,1)) &
1264  rhodotthermalannihilation(s,9) = max(- 4.0_preal * rhodip(s,1) * vclimb(s) / (dupper(s,1) - dlower(s,1)), &
1265  - rhodip(s,1) / timestep - rhodotathermalannihilation(s,9) &
1266  - rhodotsingle2dipoleglide(s,9)) ! make sure that we do not annihilate more dipoles than we have
1267 
1268  rhodot = rhodotflux &
1269  + rhodotmultiplication &
1270  + rhodotsingle2dipoleglide &
1271  + rhodotathermalannihilation &
1272  + rhodotthermalannihilation
1273 
1274 # 1325 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive_plastic_nonlocal.f90"
1275 
1276 
1277  if ( any(rho(:,mob) + rhodot(:,1:4) * timestep < -prm%atol_rho) &
1278  .or. any(rho(:,dip) + rhodot(:,9:10) * timestep < -prm%atol_rho)) then
1279 
1280 
1281 
1282 
1283 
1284 
1285  plasticstate(ph)%dotState = ieee_value(1.0_preal,ieee_quiet_nan)
1286  else
1287  dot%rho(:,of) = pack(rhodot,.true.)
1288  dot%gamma(:,of) = sum(gdot,2)
1289  endif
1290 
1291  end associate
1292 
1293 end subroutine plastic_nonlocal_dotstate
1294 
1295 
1296 !--------------------------------------------------------------------------------------------------
1299 ! plane normals and signed cosine of the angle between the slip directions. Only the largest values
1300 ! that sum up to a total of 1 are considered, all others are set to zero.
1301 !--------------------------------------------------------------------------------------------------
1302 module subroutine plastic_nonlocal_updatecompatibility(orientation,instance,i,e)
1303 
1304  type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: &
1305  orientation ! crystal orientation
1306  integer, intent(in) :: &
1307  instance, &
1308  i, &
1309  e
1310 
1311  integer :: &
1312  n, & ! neighbor index
1313  neighbor_e, & ! element index of my neighbor
1314  neighbor_i, & ! integration point index of my neighbor
1315  ph, &
1316  neighbor_phase, &
1317  ns, & ! number of active slip systems
1318  s1, & ! slip system index (me)
1319  s2 ! slip system index (my neighbor)
1320  real(pReal), dimension(2,param(instance)%sum_N_sl,param(instance)%sum_N_sl,nIPneighbors) :: &
1321  my_compatibility ! my_compatibility for current element and ip
1322  real(pReal) :: &
1323  my_compatibilitySum, &
1324  thresholdValue, &
1325  nThresholdValues
1326  logical, dimension(param(instance)%sum_N_sl) :: &
1327  belowThreshold
1328  type(rotation) :: mis
1329 
1330  ph = material_phaseat(1,e)
1331 
1332  associate(prm => param(instance))
1333  ns = prm%sum_N_sl
1334 
1335  !*** start out fully compatible
1336  my_compatibility = 0.0_preal
1337  forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_preal
1338 
1339  neighbors: do n = 1,nipneighbors
1340  neighbor_e = ipneighborhood(1,n,i,e)
1341  neighbor_i = ipneighborhood(2,n,i,e)
1342 
1343  neighbor_phase = material_phaseat(1,neighbor_e)
1344 
1345  if (neighbor_e <= 0 .or. neighbor_i <= 0) then
1346  !* FREE SURFACE
1347  !* Set surface transmissivity to the value specified in the material.config
1348  forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%surfaceTransmissivity)
1349  elseif (neighbor_phase /= ph) then
1350  !* PHASE BOUNDARY
1351  !* If we encounter a different nonlocal phase at the neighbor,
1352  !* we consider this to be a real "physical" phase boundary, so completely incompatible.
1353  !* If one of the two phases has a local plasticity law,
1354  !* we do not consider this to be a phase boundary, so completely compatible.
1355  if (.not. phase_localplasticity(neighbor_phase) .and. .not. phase_localplasticity(ph)) &
1356  forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_preal
1357  elseif (prm%grainboundaryTransmissivity >= 0.0_preal) then
1358  !* GRAIN BOUNDARY !
1359  !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
1360  if (material_texture(1,i,e) /= material_texture(1,neighbor_i,neighbor_e) .and. &
1361  (.not. phase_localplasticity(neighbor_phase))) &
1362  forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%grainboundaryTransmissivity)
1363  else
1364  !* GRAIN BOUNDARY ?
1365  !* Compatibility defined by relative orientation of slip systems:
1366  !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection.
1367  !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection.
1368  !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one),
1369  !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that
1370  !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
1371  !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
1372  !* All values below the threshold are set to zero.
1373  mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e))
1374  myslipsystems: do s1 = 1,ns
1375  neighborslipsystems: do s2 = 1,ns
1376  my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), &
1377  mis%rotate(prm%slip_normal(1:3,s2))) &
1378  * abs(math_inner(prm%slip_direction(1:3,s1), &
1379  mis%rotate(prm%slip_direction(1:3,s2))))
1380  my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
1381  mis%rotate(prm%slip_normal(1:3,s2)))) &
1382  * abs(math_inner(prm%slip_direction(1:3,s1), &
1383  mis%rotate(prm%slip_direction(1:3,s2))))
1384  enddo neighborslipsystems
1385 
1386  my_compatibilitysum = 0.0_preal
1387  belowthreshold = .true.
1388  do while (my_compatibilitysum < 1.0_preal .and. any(belowthreshold))
1389  thresholdvalue = maxval(my_compatibility(2,:,s1,n), belowthreshold) ! screws always positive
1390  nthresholdvalues = real(count(my_compatibility(2,:,s1,n) >= thresholdvalue),preal)
1391  where (my_compatibility(2,:,s1,n) >= thresholdvalue) belowthreshold = .false.
1392  if (my_compatibilitysum + thresholdvalue * nthresholdvalues > 1.0_preal) &
1393  where (abs(my_compatibility(:,:,s1,n)) >= thresholdvalue) &
1394  my_compatibility(:,:,s1,n) = sign((1.0_preal - my_compatibilitysum)/nthresholdvalues,&
1395  my_compatibility(:,:,s1,n))
1396  my_compatibilitysum = my_compatibilitysum + nthresholdvalues * thresholdvalue
1397  enddo
1398 
1399  where(belowthreshold) my_compatibility(1,:,s1,n) = 0.0_preal
1400  where(belowthreshold) my_compatibility(2,:,s1,n) = 0.0_preal
1401 
1402  enddo myslipsystems
1403  endif
1404 
1405  enddo neighbors
1406 
1407  compatibility(:,:,:,:,i,e) = my_compatibility
1408 
1409  end associate
1410 
1411 end subroutine plastic_nonlocal_updatecompatibility
1412 
1413 
1414 !--------------------------------------------------------------------------------------------------
1416 !--------------------------------------------------------------------------------------------------
1417 module subroutine plastic_nonlocal_results(instance,group)
1418 
1419  integer, intent(in) :: instance
1420  character(len=*),intent(in) :: group
1421 
1422  integer :: o
1423 
1424  associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance))
1425  outputsloop: do o = 1,size(prm%output)
1426  select case(trim(prm%output(o)))
1427  case('rho_sgl_mob_edg_pos')
1428  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', &
1429  'positive mobile edge density','1/m²')
1430  case('rho_sgl_imm_edg_pos')
1431  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',&
1432  'positive immobile edge density','1/m²')
1433  case('rho_sgl_mob_edg_neg')
1434  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',&
1435  'negative mobile edge density','1/m²')
1436  case('rho_sgl_imm_edg_neg')
1437  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',&
1438  'negative immobile edge density','1/m²')
1439  case('rho_dip_edg')
1440  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_dip_edg, 'rho_dip_edg',&
1441  'edge dipole density','1/m²')
1442  case('rho_sgl_mob_scr_pos')
1443  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',&
1444  'positive mobile screw density','1/m²')
1445  case('rho_sgl_imm_scr_pos')
1446  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',&
1447  'positive immobile screw density','1/m²')
1448  case('rho_sgl_mob_scr_neg')
1449  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',&
1450  'negative mobile screw density','1/m²')
1451  case('rho_sgl_imm_scr_neg')
1452  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',&
1453  'negative immobile screw density','1/m²')
1454  case('rho_dip_scr')
1455  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_dip_scr, 'rho_dip_scr',&
1456  'screw dipole density','1/m²')
1457  case('rho_forest')
1458  if(prm%sum_N_sl>0) call results_writedataset(group,stt%rho_forest, 'rho_forest',&
1459  'forest density','1/m²')
1460  case('v_edg_pos')
1461  if(prm%sum_N_sl>0) call results_writedataset(group,stt%v_edg_pos, 'v_edg_pos',&
1462  'positive edge velocity','m/s')
1463  case('v_edg_neg')
1464  if(prm%sum_N_sl>0) call results_writedataset(group,stt%v_edg_neg, 'v_edg_neg',&
1465  'negative edge velocity','m/s')
1466  case('v_scr_pos')
1467  if(prm%sum_N_sl>0) call results_writedataset(group,stt%v_scr_pos, 'v_scr_pos',&
1468  'positive srew velocity','m/s')
1469  case('v_scr_neg')
1470  if(prm%sum_N_sl>0) call results_writedataset(group,stt%v_scr_neg, 'v_scr_neg',&
1471  'negative screw velocity','m/s')
1472  case('gamma')
1473  if(prm%sum_N_sl>0) call results_writedataset(group,stt%gamma,'gamma',&
1474  'plastic shear','1')
1475  case('tau_pass')
1476  if(prm%sum_N_sl>0) call results_writedataset(group,dst%tau_pass,'tau_pass',&
1477  'passing stress for slip','Pa')
1478  end select
1479  enddo outputsloop
1480  end associate
1481 
1482 end subroutine plastic_nonlocal_results
1483 
1484 
1485 !--------------------------------------------------------------------------------------------------
1487 !--------------------------------------------------------------------------------------------------
1488 subroutine stateinit(ini,phase,NipcMyPhase)
1489 
1490  type(tInitialParameters) :: &
1491  ini
1492  integer,intent(in) :: &
1493  phase, &
1494  NipcMyPhase
1495  integer :: &
1496  e, &
1497  i, &
1498  f, &
1499  from, &
1500  upto, &
1501  s, &
1502  instance, &
1503  phasemember
1504  real(pReal), dimension(2) :: &
1505  noise, &
1506  rnd
1507  real(pReal) :: &
1508  meanDensity, &
1509  totalVolume, &
1510  densityBinning, &
1511  minimumIpVolume
1512  real(pReal), dimension(NipcMyPhase) :: &
1513  volume
1514 
1515  instance = phase_plasticityinstance(phase)
1516  associate(stt => state(instance))
1517 
1518  if (ini%rhoSglRandom > 0.0_preal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
1519  do e = 1,discretization_nelem
1520  do i = 1,discretization_nip
1521  if (material_phaseat(1,e) == phase) volume(material_phasememberat(1,i,e)) = ipvolume(i,e)
1522  enddo
1523  enddo
1524  totalvolume = sum(volume)
1525  minimumipvolume = minval(volume)
1526  densitybinning = ini%rhoSglRandomBinning / minimumipvolume ** (2.0_preal / 3.0_preal)
1527 
1528  ! fill random material points with dislocation segments until the desired overall density is reached
1529  meandensity = 0.0_preal
1530  do while(meandensity < ini%rhoSglRandom)
1531  call random_number(rnd)
1532  phasemember = nint(rnd(1)*real(nipcmyphase,preal) + 0.5_preal)
1533  s = nint(rnd(2)*real(sum(ini%N_sl),preal)*4.0_preal + 0.5_preal)
1534  meandensity = meandensity + densitybinning * volume(phasemember) / totalvolume
1535  stt%rhoSglMobile(s,phasemember) = densitybinning
1536  enddo
1537  else ! homogeneous distribution with noise
1538  do e = 1, nipcmyphase
1539  do f = 1,size(ini%N_sl,1)
1540  from = 1 + sum(ini%N_sl(1:f-1))
1541  upto = sum(ini%N_sl(1:f))
1542  do s = from,upto
1543  noise = [math_samplegaussvar(0.0_preal, ini%rhoSglScatter), &
1544  math_samplegaussvar(0.0_preal, ini%rhoSglScatter)]
1545  stt%rho_sgl_mob_edg_pos(s,e) = ini%rhoSglEdgePos0(f) + noise(1)
1546  stt%rho_sgl_mob_edg_neg(s,e) = ini%rhoSglEdgeNeg0(f) + noise(1)
1547  stt%rho_sgl_mob_scr_pos(s,e) = ini%rhoSglScrewPos0(f) + noise(2)
1548  stt%rho_sgl_mob_scr_neg(s,e) = ini%rhoSglScrewNeg0(f) + noise(2)
1549  enddo
1550  stt%rho_dip_edg(from:upto,e) = ini%rhoDipEdge0(f)
1551  stt%rho_dip_scr(from:upto,e) = ini%rhoDipScrew0(f)
1552  enddo
1553  enddo
1554  endif
1555 
1556  end associate
1557 
1558 end subroutine stateinit
1559 
1560 
1561 !--------------------------------------------------------------------------------------------------
1563 !--------------------------------------------------------------------------------------------------
1564 subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance)
1565 
1566  integer, intent(in) :: &
1567  c, & !< dislocation character (1:edge, 2:screw)
1568  instance
1569  real(pReal), intent(in) :: &
1570  Temperature
1571  real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: &
1572  tau, & !< resolved external shear stress (without non Schmid effects)
1573  tauNS, & !< resolved external shear stress (including non Schmid effects)
1574  tauThreshold
1575  real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: &
1576  v, & !< velocity
1577  dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions)
1578  dv_dtauNS
1579 
1580  integer :: &
1581  ns, & !< short notation for the total number of active slip systems
1582  s
1583  real(pReal) :: &
1584  tauRel_P, &
1585  tauRel_S, &
1586  tauEff, & !< effective shear stress
1587  tPeierls, & !< waiting time in front of a peierls barriers
1588  tSolidSolution, & !< waiting time in front of a solid solution obstacle
1589  vViscous, & !< viscous glide velocity
1590  dtPeierls_dtau, & !< derivative with respect to resolved shear stress
1591  dtSolidSolution_dtau, & !< derivative with respect to resolved shear stress
1592  meanfreepath_S, & !< mean free travel distance for dislocations between two solid solution obstacles
1593  meanfreepath_P, & !< mean free travel distance for dislocations between two Peierls barriers
1594  jumpWidth_P, & !< depth of activated area
1595  jumpWidth_S, & !< depth of activated area
1596  activationLength_P, & !< length of activated dislocation line
1597  activationLength_S, & !< length of activated dislocation line
1598  activationVolume_P, & !< volume that needs to be activated to overcome barrier
1599  activationVolume_S, & !< volume that needs to be activated to overcome barrier
1600  activationEnergy_P, & !< energy that is needed to overcome barrier
1601  activationEnergy_S, & !< energy that is needed to overcome barrier
1602  criticalStress_P, & !< maximum obstacle strength
1603  criticalStress_S, & !< maximum obstacle strength
1604  mobility
1605 
1606  associate(prm => param(instance))
1607  ns = prm%sum_N_sl
1608  v = 0.0_preal
1609  dv_dtau = 0.0_preal
1610  dv_dtauns = 0.0_preal
1611 
1612  do s = 1,ns
1613  if (abs(tau(s)) > tauthreshold(s)) then
1614 
1615  !* Peierls contribution
1616  !* Effective stress includes non Schmid constributions
1617  !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
1618  taueff = max(0.0_preal, abs(tauns(s)) - tauthreshold(s)) ! ensure that the effective stress is positive
1619  meanfreepath_p = prm%burgers(s)
1620  jumpwidth_p = prm%burgers(s)
1621  activationlength_p = prm%doublekinkwidth *prm%burgers(s)
1622  activationvolume_p = activationlength_p * jumpwidth_p * prm%burgers(s)
1623  criticalstress_p = prm%peierlsStress(s,c)
1624  activationenergy_p = criticalstress_p * activationvolume_p
1625  taurel_p = min(1.0_preal, taueff / criticalstress_p) ! ensure that the activation probability cannot become greater than one
1626  tpeierls = 1.0_preal / prm%fattack &
1627  * exp(activationenergy_p / (kb * temperature) &
1628  * (1.0_preal - taurel_p**prm%p)**prm%q)
1629  if (taueff < criticalstress_p) then
1630  dtpeierls_dtau = tpeierls * prm%p * prm%q * activationvolume_p / (kb * temperature) &
1631  * (1.0_preal - taurel_p**prm%p)**(prm%q-1.0_preal) * taurel_p**(prm%p-1.0_preal)
1632  else
1633  dtpeierls_dtau = 0.0_preal
1634  endif
1635 
1636  !* Contribution from solid solution strengthening
1637  !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
1638  taueff = abs(tau(s)) - tauthreshold(s)
1639  meanfreepath_s = prm%burgers(s) / sqrt(prm%solidSolutionConcentration)
1640  jumpwidth_s = prm%solidSolutionSize * prm%burgers(s)
1641  activationlength_s = prm%burgers(s) / sqrt(prm%solidSolutionConcentration)
1642  activationvolume_s = activationlength_s * jumpwidth_s * prm%burgers(s)
1643  activationenergy_s = prm%solidSolutionEnergy
1644  criticalstress_s = activationenergy_s / activationvolume_s
1645  taurel_s = min(1.0_preal, taueff / criticalstress_s) ! ensure that the activation probability cannot become greater than one
1646  tsolidsolution = 1.0_preal / prm%fattack &
1647  * exp(activationenergy_s / (kb * temperature)* (1.0_preal - taurel_s**prm%p)**prm%q)
1648  if (taueff < criticalstress_s) then
1649  dtsolidsolution_dtau = tsolidsolution * prm%p * prm%q * activationvolume_s / (kb * temperature) &
1650  * (1.0_preal - taurel_s**prm%p)**(prm%q-1.0_preal)* taurel_s**(prm%p-1.0_preal)
1651  else
1652  dtsolidsolution_dtau = 0.0_preal
1653  endif
1654 
1655  !* viscous glide velocity
1656  taueff = abs(tau(s)) - tauthreshold(s)
1657  mobility = prm%burgers(s) / prm%viscosity
1658  vviscous = mobility * taueff
1659 
1660  !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of
1661  !* free flight at glide velocity in between.
1662  !* adopt sign from resolved stress
1663  v(s) = sign(1.0_preal,tau(s)) &
1664  / (tpeierls / meanfreepath_p + tsolidsolution / meanfreepath_s + 1.0_preal / vviscous)
1665  dv_dtau(s) = v(s)**2.0_preal * (dtsolidsolution_dtau / meanfreepath_s + mobility /vviscous**2.0_preal)
1666  dv_dtauns(s) = v(s)**2.0_preal * dtpeierls_dtau / meanfreepath_p
1667  endif
1668  enddo
1669 
1670  end associate
1671 
1672 end subroutine kinetics
1673 
1674 
1675 !--------------------------------------------------------------------------------------------------
1678 !--------------------------------------------------------------------------------------------------
1679 function getrho(instance,of,ip,el)
1680 
1681  integer, intent(in) :: instance, of,ip,el
1682  real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho
1683 
1684  associate(prm => param(instance))
1685 
1686  getrho = reshape(state(instance)%rho(:,of),[prm%sum_N_sl,10])
1687 
1688  ! ensure positive densities (not for imm, they have a sign)
1689  getrho(:,mob) = max(getrho(:,mob),0.0_preal)
1690  getrho(:,dip) = max(getrho(:,dip),0.0_preal)
1691 
1692  where(abs(getrho) < max(prm%significantN/ipvolume(ip,el)**(2.0_preal/3.0_preal),prm%significantRho)) &
1693  getrho = 0.0_preal
1694 
1695  end associate
1696 
1697 end function getrho
1698 
1699 
1700 !--------------------------------------------------------------------------------------------------
1703 !--------------------------------------------------------------------------------------------------
1704 function getrho0(instance,of,ip,el)
1705 
1706  integer, intent(in) :: instance, of,ip,el
1707  real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho0
1708 
1709  associate(prm => param(instance))
1710 
1711  getrho0 = reshape(state0(instance)%rho(:,of),[prm%sum_N_sl,10])
1712 
1713  ! ensure positive densities (not for imm, they have a sign)
1714  getrho0(:,mob) = max(getrho0(:,mob),0.0_preal)
1715  getrho0(:,dip) = max(getrho0(:,dip),0.0_preal)
1716 
1717  where(abs(getrho0) < max(prm%significantN/ipvolume(ip,el)**(2.0_preal/3.0_preal),prm%significantRho)) &
1718  getrho0 = 0.0_preal
1719 
1720  end associate
1721 
1722 end function getrho0
1723 
1724 end submodule plastic_nonlocal
geometry_plastic_nonlocal::geometry_plastic_nonlocal_ipareanormal0
real(preal), dimension(:,:,:,:), allocatable, protected geometry_plastic_nonlocal_ipareanormal0
area normal of interface to neighboring IP (initially!)
Definition: geometry_plastic_nonlocal.f90:31
geometry_plastic_nonlocal::geometry_plastic_nonlocal_ipvolume0
real(preal), dimension(:,:), allocatable, protected geometry_plastic_nonlocal_ipvolume0
volume associated with IP (initially!)
Definition: geometry_plastic_nonlocal.f90:25
geometry_plastic_nonlocal::geometry_plastic_nonlocal_nipneighbors
integer, protected geometry_plastic_nonlocal_nipneighbors
Definition: geometry_plastic_nonlocal.f90:19
config
Reads in the material configuration from file.
Definition: config.f90:13
geometry_plastic_nonlocal
Geometric information about the IP cells needed for the nonlocal.
Definition: geometry_plastic_nonlocal.f90:12
geometry_plastic_nonlocal::geometry_plastic_nonlocal_ipneighborhood
integer, dimension(:,:,:,:), allocatable, protected geometry_plastic_nonlocal_ipneighborhood
6 or less neighboring IPs as [element ID, IP ID, face ID that point to me]
Definition: geometry_plastic_nonlocal.f90:22
geometry_plastic_nonlocal::geometry_plastic_nonlocal_iparea0
real(preal), dimension(:,:,:), allocatable, protected geometry_plastic_nonlocal_iparea0
area of interface to neighboring IP (initially!)
Definition: geometry_plastic_nonlocal.f90:28
constitutive
elasticity, plasticity, internal microstructure state
Definition: constitutive.f90:10