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

Merge branch 'development' into misc-improvements

This commit is contained in:
Martin Diehl 2020-04-24 20:01:40 +02:00
parent cc3fa156cf 05fd943cf4
commit 1610d5a5d2
36 changed files with 531 additions and 588 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitlab-ci.yml
View File

@ -203,7 +203,6 @@ Post_OrientationConversion:
stage: postprocessing
script:
- OrientationConversion/test.py
- OrientationConversion/test2.py
except:
- master
- release

2
VERSION
View File

@ -1 +1 @@
v2.0.3-2303-g2a6132b7
v2.0.3-2364-g62f7363a

2
processing/post/DADF5_postResults.py
View File

@ -33,7 +33,7 @@ for filename in options.filenames:
results = damask.Result(filename)
if not results.structured: continue
coords = damask.grid_filters.cell_coord0(results.grid,results.size,results.origin)
coords = damask.grid_filters.cell_coord0(results.grid,results.size,results.origin).reshape(-1,3,order='F')
N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1
N_digits = 5 # hack to keep test intact

55
processing/post/addCompatibilityMismatch.py
View File

@ -17,7 +17,7 @@ def volTetrahedron(coords):
"""
Return the volume of the tetrahedron with given vertices or sides.
Ifvertices are given they must be in a NumPy array with shape (4,3): the
If vertices are given they must be in a NumPy array with shape (4,3): the
position vectors of the 4 vertices in 3 dimensions; if the six sides are
given, they must be an array of length 6. If both are given, the sides
will be used in the calculation.
@ -67,14 +67,13 @@ def volumeMismatch(size,F,nodes):
(compatible) cube and determinant of deformation gradient at Fourier point.
"""
coords = np.empty([8,3])
vMismatch = np.empty(grid[::-1])
volInitial = size.prod()/grid.prod()
vMismatch = np.empty(F.shape[:3])
#--------------------------------------------------------------------------------------------------
# calculate actual volume and volume resulting from deformation gradient
for k in range(grid[2]):
for k in range(grid[0]):
for j in range(grid[1]):
for i in range(grid[0]):
for i in range(grid[2]):
coords[0,0:3] = nodes[k, j, i ,0:3]
coords[1,0:3] = nodes[k ,j, i+1,0:3]
coords[2,0:3] = nodes[k ,j+1,i+1,0:3]
@ -91,8 +90,7 @@ def volumeMismatch(size,F,nodes):
+ abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[5,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[0,0:3]]))) \
/np.linalg.det(F[k,j,i,0:3,0:3])
return vMismatch/volInitial
return vMismatch/(size.prod()/grid.prod())
def shapeMismatch(size,F,nodes,centres):
@ -103,35 +101,34 @@ def shapeMismatch(size,F,nodes,centres):
the corners of reconstructed (combatible) volume element and the vectors calculated by deforming
the initial volume element with the current deformation gradient.
"""
coordsInitial = np.empty([8,3])
sMismatch = np.empty(grid[::-1])
sMismatch = np.empty(F.shape[:3])
#--------------------------------------------------------------------------------------------------
# initial positions
coordsInitial[0,0:3] = [-size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[1,0:3] = [+size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[2,0:3] = [+size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[3,0:3] = [-size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[4,0:3] = [-size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
coordsInitial[5,0:3] = [+size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
coordsInitial[6,0:3] = [+size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
coordsInitial[7,0:3] = [-size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
coordsInitial = coordsInitial/2.0
delta = size/grid*.5
coordsInitial = np.vstack((delta * np.array((-1,-1,-1)),
delta * np.array((+1,-1,-1)),
delta * np.array((+1,+1,-1)),
delta * np.array((-1,+1,-1)),
delta * np.array((-1,-1,+1)),
delta * np.array((+1,-1,+1)),
delta * np.array((+1,+1,+1)),
delta * np.array((-1,+1,+1))))
#--------------------------------------------------------------------------------------------------
# compare deformed original and deformed positions to actual positions
for k in range(grid[2]):
for k in range(grid[0]):
for j in range(grid[1]):
for i in range(grid[0]):
for i in range(grid[2]):
sMismatch[k,j,i] = \
+ np.linalg.norm(nodes[k, j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[0,0:3]))\
+ np.linalg.norm(nodes[k, j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\
+ np.linalg.norm(nodes[k, j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\
+ np.linalg.norm(nodes[k+1,j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\
+ np.linalg.norm(nodes[k+1,j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\
+ np.linalg.norm(nodes[k, j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[3,0:3]))\
+ np.linalg.norm(nodes[k+1,j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\
+ np.linalg.norm(nodes[k, j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\
+ np.linalg.norm(nodes[k+1,j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[5,0:3]))\
+ np.linalg.norm(nodes[k+1,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[6,0:3]))\
+ np.linalg.norm(nodes[k+1,j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3]))
+ np.linalg.norm(nodes[k ,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3]))
return sMismatch
@ -178,20 +175,20 @@ for name in filenames:
table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
F = table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3)
F = table.get(options.defgrad).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3))
nodes = damask.grid_filters.node_coord(size,F)
if options.shape:
centers = damask.grid_filters.cell_coord(size,F)
shapeMismatch = shapeMismatch( size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes,centers)
shapeMismatch = shapeMismatch(size,F,nodes,centers)
table.add('shapeMismatch(({}))'.format(options.defgrad),
shapeMismatch.reshape(-1,1),
shapeMismatch.reshape(-1,1,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))
if options.volume:
volumeMismatch = volumeMismatch(size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes)
volumeMismatch = volumeMismatch(size,F,nodes)
table.add('volMismatch(({}))'.format(options.defgrad),
volumeMismatch.reshape(-1,1),
volumeMismatch.reshape(-1,1,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name)

5
processing/post/addCurl.py
View File

@ -49,9 +49,10 @@ for name in filenames:
for label in options.labels:
field = table.get(label)
shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor
field = field.reshape(np.append(grid[::-1],shape))
field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape)
curl = damask.grid_filters.curl(size,field)
table.add('curlFFT({})'.format(label),
damask.grid_filters.curl(size[::-1],field).reshape(-1,np.prod(shape)),
curl.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape),order='F'),
scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name)

14
processing/post/addDisplacement.py
View File

@ -5,8 +5,6 @@ import sys
from io import StringIO
from optparse import OptionParser
import numpy as np
import damask
@ -52,22 +50,22 @@ for name in filenames:
table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
F = table.get(options.f).reshape(np.append(grid[::-1],(3,3)))
F = table.get(options.f).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3))
if options.nodal:
table = damask.Table(damask.grid_filters.node_coord0(grid[::-1],size[::-1]).reshape(-1,3),
table = damask.Table(damask.grid_filters.node_coord0(grid,size).reshape(-1,3,order='F'),
{'pos':(3,)})
table.add('avg({}).{}'.format(options.f,options.pos),
damask.grid_filters.node_displacement_avg(size[::-1],F).reshape(-1,3),
damask.grid_filters.node_displacement_avg(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))
table.add('fluct({}).{}'.format(options.f,options.pos),
damask.grid_filters.node_displacement_fluct(size[::-1],F).reshape(-1,3),
damask.grid_filters.node_displacement_fluct(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt')
else:
table.add('avg({}).{}'.format(options.f,options.pos),
damask.grid_filters.cell_displacement_avg(size[::-1],F).reshape(-1,3),
damask.grid_filters.cell_displacement_avg(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))
table.add('fluct({}).{}'.format(options.f,options.pos),
damask.grid_filters.cell_displacement_fluct(size[::-1],F).reshape(-1,3),
damask.grid_filters.cell_displacement_fluct(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name)

5
processing/post/addDivergence.py
View File

@ -49,9 +49,10 @@ for name in filenames:
for label in options.labels:
field = table.get(label)
shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor
field = field.reshape(np.append(grid[::-1],shape))
field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape)
div = damask.grid_filters.divergence(size,field)
table.add('divFFT({})'.format(label),
damask.grid_filters.divergence(size[::-1],field).reshape(-1,np.prod(shape)//3),
div.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)//3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name)

5
processing/post/addGradient.py
View File

@ -49,9 +49,10 @@ for name in filenames:
for label in options.labels:
field = table.get(label)
shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,) # scalar or vector
field = field.reshape(np.append(grid[::-1],shape))
field = field.reshape(tuple(grid)+(-1,),order='F')
grad = damask.grid_filters.gradient(size,field)
table.add('gradFFT({})'.format(label),
damask.grid_filters.gradient(size[::-1],field).reshape(-1,np.prod(shape)*3),
grad.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)*3,order='F'),
scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name)

2
processing/post/averageDown.py
View File

@ -91,7 +91,7 @@ for name in filenames:
table = damask.Table(averagedDown,table.shapes,table.comments)
coords = damask.grid_filters.cell_coord0(packedGrid,size,shift/packedGrid*size+origin)
table.set(options.pos, coords.reshape(-1,3))
table.set(options.pos, coords.reshape(-1,3,order='F'))
outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name))

6
processing/post/blowUp.py
View File

@ -59,13 +59,13 @@ for name in filenames:
packing = np.array(options.packing,'i')
outSize = grid*packing
data = table.data.values.reshape(tuple(grid)+(-1,))
blownUp = ndimage.interpolation.zoom(data,tuple(packing)+(1,),order=0,mode='nearest').reshape(outSize.prod(),-1)
data = table.data.values.reshape(tuple(grid)+(-1,),order='F')
blownUp = ndimage.interpolation.zoom(data,tuple(packing)+(1,),order=0,mode='nearest').reshape(outSize.prod(),-1,order='F')
table = damask.Table(blownUp,table.shapes,table.comments)
coords = damask.grid_filters.cell_coord0(outSize,size,origin)
table.set(options.pos,coords.reshape(-1,3))
table.set(options.pos,coords.reshape(-1,3,order='F'))
table.set('elem',np.arange(1,outSize.prod()+1))
outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name))

6
processing/post/permuteData.py
View File

@ -97,10 +97,10 @@ for name in filenames:
for col,dim in zip(columns,dims):
if options.unique:
s = set(map(tuple,table.data[:,col:col+dim])) # generate set of (unique) values
uniques = np.array(map(np.array,s)) # translate set to np.array
uniques = np.array(list(map(np.array,s))) # translate set to np.array
shuffler = dict(zip(s,np.random.permutation(len(s)))) # random permutation
table.data[:,col:col+dim] = uniques[np.array(map(lambda x: shuffler[tuple(x)],
table.data[:,col:col+dim]))] # fill table with mapped uniques
table.data[:,col:col+dim] = uniques[np.array(list(map(lambda x: shuffler[tuple(x)],
table.data[:,col:col+dim])))] # fill table with mapped uniques
else:
np.random.shuffle(table.data[:,col:col+dim]) # independently shuffle every row

12
processing/pre/geom_fromVoronoiTessellation.py
View File

@ -24,22 +24,22 @@ def findClosestSeed(seeds, weights, point):
def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), periodic = True, cpus = 2):
if periodic:
weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
weights_p = np.tile(weights.squeeze(),27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3,order='F')
coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3)
else:
weights_p = weights.flatten()
weights_p = weights.squeeze()
seeds_p = seeds
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3)
if cpus > 1:
pool = multiprocessing.Pool(processes = cpus)
result = pool.map_async(partial(findClosestSeed,seeds_p,weights_p), [coord for coord in coords])
pool.close()
pool.join()
closest_seed = np.array(result.get())
closest_seed = np.array(result.get()).reshape(-1,3)
else:
closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords])
@ -52,7 +52,7 @@ def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), peri
def Voronoi_tessellation(grid, size, seeds, origin = np.zeros(3), periodic = True):
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,closest_seed = KDTree.query(coords)

2
processing/pre/seeds_fromGeom.py
View File

@ -54,7 +54,7 @@ for name in filenames:
np.in1d(microstructure,options.blacklist,invert=True) if options.blacklist else \
np.full(geom.grid.prod(),True,dtype=bool))
seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3)
seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
comments = geom.comments \
+ [scriptID + ' ' + ' '.join(sys.argv[1:]),

2
processing/pre/seeds_fromRandom.py
View File

@ -128,7 +128,7 @@ for name in filenames:
if not options.selective:
coords = damask.grid_filters.cell_coord0(grid,size).reshape(-1,3)
coords = damask.grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
seeds = coords[np.random.choice(np.prod(grid), options.N, replace=False)] \
+ np.broadcast_to(size/grid,(options.N,3))*(np.random.rand(options.N,3)*.5-.25) # wobble without leaving grid
else:

16
python/damask/_geom.py
View File

@ -322,11 +322,10 @@ class Geom:
if i != grid.prod():
raise TypeError('Invalid file: expected {} entries, found {}'.format(grid.prod(),i))
microstructure = microstructure.reshape(grid,order='F')
if not np.any(np.mod(microstructure.flatten(),1) != 0.0): # no float present
if not np.any(np.mod(microstructure,1) != 0.0): # no float present
microstructure = microstructure.astype('int')
return Geom(microstructure.reshape(grid),size,origin,homogenization,comments)
return Geom(microstructure.reshape(grid,order='F'),size,origin,homogenization,comments)
@staticmethod
@ -352,16 +351,15 @@ class Geom:
"""
if periodic:
weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F')
coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3)
else:
weights_p = weights.flatten()
weights_p = weights
seeds_p = seeds
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
pool = multiprocessing.Pool(processes = int(Environment().options['DAMASK_NUM_THREADS']))
result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
@ -396,7 +394,7 @@ class Geom:
perform a periodic tessellation. Defaults to True.
"""
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,microstructure = KDTree.query(coords)

16
python/damask/_result.py
View File

@ -111,7 +111,7 @@ class Result:
select from 'set', 'add', and 'del'
what : str
attribute to change (must be from self.selection)
datasets : list of str or Boolean
datasets : list of str or bool
name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to []
@ -203,7 +203,7 @@ class Result:
----------
what : str
attribute to change (must be from self.selection)
datasets : list of str or Boolean
datasets : list of str or bool
name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to []
@ -219,7 +219,7 @@ class Result:
----------
what : str
attribute to change (must be from self.selection)
datasets : list of str or Boolean
datasets : list of str or bool
name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to []
@ -235,7 +235,7 @@ class Result:
----------
what : str
attribute to change (must be from self.selection)
datasets : list of str or Boolean
datasets : list of str or bool
name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to []
@ -262,10 +262,10 @@ class Result:
datasets : iterable or str
component : int
homogenization component to consider for constituent data
tagged : Boolean
tagged : bool
tag Table.column name with '#component'
defaults to False
split : Boolean
split : bool
split Table by increment and return dictionary of Tables
defaults to True
@ -326,7 +326,7 @@ class Result:
Parameters
----------
datasets : iterable or str or Boolean
datasets : iterable or str or bool
Examples
--------
@ -460,7 +460,7 @@ class Result:
def cell_coordinates(self):
"""Return initial coordinates of the cell centers."""
if self.structured:
return grid_filters.cell_coord0(self.grid,self.size,self.origin).reshape(-1,3)
return grid_filters.cell_coord0(self.grid,self.size,self.origin).reshape(-1,3,order='F')
else:
with h5py.File(self.fname,'r') as f:
return f['geometry/x_c'][()]

132
python/damask/grid_filters.py
View File

@ -1,3 +1,17 @@
"""
Filters for operations on regular grids.
Notes
-----
The grids are defined as (x,y,z,...) where x is fastest and z is slowest.
This convention is consistent with the geom file format.
When converting to/from a plain list (e.g. storage in ASCII table),
the following operations are required for tensorial data:
D3 = D1.reshape(grid+(-1,),order='F').reshape(grid+(3,3))
D1 = D3.reshape(grid+(-1,)).reshape(-1,9,order='F')
"""
from scipy import spatial as _spatial
import numpy as _np
@ -7,8 +21,12 @@ def _ks(size,grid,first_order=False):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
grid : numpy.ndarray of shape (3)
number of grid points.
first_order : bool, optional
correction for first order derivatives, defaults to False.
"""
k_sk = _np.where(_np.arange(grid[0])>grid[0]//2,_np.arange(grid[0])-grid[0],_np.arange(grid[0]))/size[0]
@ -19,8 +37,7 @@ def _ks(size,grid,first_order=False):
k_si = _np.arange(grid[2]//2+1)/size[2]
kk, kj, ki = _np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
return _np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3)
return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
def curl(size,field):
@ -29,8 +46,10 @@ def curl(size,field):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
periodic field of which the curl is calculated.
"""
n = _np.prod(field.shape[3:])
@ -53,8 +72,10 @@ def divergence(size,field):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
periodic field of which the divergence is calculated.
"""
n = _np.prod(field.shape[3:])
@ -69,12 +90,14 @@ def divergence(size,field):
def gradient(size,field):
"""
Calculate gradient of a vector or scalar field in Fourier space.
Calculate gradient of a scalar or vector field in Fourier space.
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3)
periodic field of which the gradient is calculated.
"""
n = _np.prod(field.shape[3:])
@ -93,9 +116,9 @@ def cell_coord0(grid,size,origin=_np.zeros(3)):
Parameters
----------
grid : numpy.ndarray
grid : numpy.ndarray of shape (3)
number of grid points.
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
origin : numpy.ndarray, optional
physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
@ -103,7 +126,11 @@ def cell_coord0(grid,size,origin=_np.zeros(3)):
"""
start = origin + size/grid*.5
end = origin + size - size/grid*.5
return _np.mgrid[start[0]:end[0]:grid[0]*1j,start[1]:end[1]:grid[1]*1j,start[2]:end[2]:grid[2]*1j].T
return _np.stack(_np.meshgrid(_np.linspace(start[0],end[0],grid[0]),
_np.linspace(start[1],end[1],grid[1]),
_np.linspace(start[2],end[2],grid[2]),indexing = 'ij'),
axis = -1)
def cell_displacement_fluct(size,F):
@ -112,7 +139,7 @@ def cell_displacement_fluct(size,F):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
F : numpy.ndarray
deformation gradient field.
@ -139,14 +166,14 @@ def cell_displacement_avg(size,F):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
F : numpy.ndarray
deformation gradient field.
"""
F_avg = _np.average(F,axis=(0,1,2))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),cell_coord0(F.shape[:3][::-1],size))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),cell_coord0(F.shape[:3],size))
def cell_displacement(size,F):
@ -155,7 +182,7 @@ def cell_displacement(size,F):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
F : numpy.ndarray
deformation gradient field.
@ -170,25 +197,25 @@ def cell_coord(size,F,origin=_np.zeros(3)):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
F : numpy.ndarray
deformation gradient field.
origin : numpy.ndarray, optional
origin : numpy.ndarray of shape (3), optional
physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
"""
return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F)
return cell_coord0(F.shape[:3],size,origin) + cell_displacement(size,F)
def cell_coord0_gridSizeOrigin(coord0,ordered=True):
"""
Return grid 'DNA', i.e. grid, size, and origin from array of cell positions.
Return grid 'DNA', i.e. grid, size, and origin from 1D array of cell positions.
Parameters
----------
coord0 : numpy.ndarray
array of undeformed cell coordinates.
coord0 : numpy.ndarray of shape (:,3)
undeformed cell coordinates.
ordered : bool, optional
expect coord0 data to be ordered (x fast, z slow).
@ -211,14 +238,14 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
start = origin + delta*.5
end = origin - delta*.5 + size
atol = 1e-4*_np.max(size)
if not _np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0]),atol=atol) and \
_np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1]),atol=atol) and \
_np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2]),atol=atol):
atol = _np.max(size)
if not (_np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0]),atol=atol) and \
_np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1]),atol=atol) and \
_np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2]),atol=atol)):
raise ValueError('Regular grid spacing violated.')
if ordered and not _np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin),atol=atol):
raise ValueError('Input data is not a regular grid.')
if ordered and not _np.allclose(coord0.reshape(tuple(grid)+(3,),order='F'),cell_coord0(grid,size,origin),atol=atol):
raise ValueError('Input data is not ordered (x fast, z slow).')
return (grid,size,origin)
@ -242,17 +269,18 @@ def node_coord0(grid,size,origin=_np.zeros(3)):
Parameters
----------
grid : numpy.ndarray
grid : numpy.ndarray of shape (3)
number of grid points.
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
origin : numpy.ndarray, optional
origin : numpy.ndarray of shape (3), optional
physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
"""
return _np.mgrid[origin[0]:size[0]+origin[0]:(grid[0]+1)*1j,
origin[1]:size[1]+origin[1]:(grid[1]+1)*1j,
origin[2]:size[2]+origin[2]:(grid[2]+1)*1j].T
return _np.stack(_np.meshgrid(_np.linspace(origin[0],size[0]+origin[0],grid[0]+1),
_np.linspace(origin[1],size[1]+origin[1],grid[1]+1),
_np.linspace(origin[2],size[2]+origin[2],grid[2]+1),indexing = 'ij'),
axis = -1)
def node_displacement_fluct(size,F):
@ -261,7 +289,7 @@ def node_displacement_fluct(size,F):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
F : numpy.ndarray
deformation gradient field.
@ -276,14 +304,14 @@ def node_displacement_avg(size,F):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
F : numpy.ndarray
deformation gradient field.
"""
F_avg = _np.average(F,axis=(0,1,2))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),node_coord0(F.shape[:3][::-1],size))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),node_coord0(F.shape[:3],size))
def node_displacement(size,F):
@ -292,7 +320,7 @@ def node_displacement(size,F):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
F : numpy.ndarray
deformation gradient field.
@ -307,15 +335,15 @@ def node_coord(size,F,origin=_np.zeros(3)):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size of the periodic field.
F : numpy.ndarray
deformation gradient field.
origin : numpy.ndarray, optional
origin : numpy.ndarray of shape (3), optional
physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
"""
return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F)
return node_coord0(F.shape[:3],size,origin) + node_displacement(size,F)
def cell_2_node(cell_data):
@ -336,14 +364,14 @@ def node_2_cell(node_data):
return c[:-1,:-1,:-1]
def node_coord0_gridSizeOrigin(coord0,ordered=False):
def node_coord0_gridSizeOrigin(coord0,ordered=True):
"""
Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions.
Return grid 'DNA', i.e. grid, size, and origin from 1D array of nodal positions.
Parameters
----------
coord0 : numpy.ndarray
array of undeformed nodal coordinates.
coord0 : numpy.ndarray of shape (:,3)
undeformed nodal coordinates.
ordered : bool, optional
expect coord0 data to be ordered (x fast, z slow).
@ -359,13 +387,13 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False):
raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid))
atol = _np.max(size)
if not _np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1),atol=atol) and \
_np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1),atol=atol) and \
_np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1),atol=atol):
if not (_np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1),atol=atol) and \
_np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1),atol=atol) and \
_np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1),atol=atol)):
raise ValueError('Regular grid spacing violated.')
if ordered and not _np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin),atol=atol):
raise ValueError('Input data is not a regular grid.')
if ordered and not _np.allclose(coord0.reshape(tuple(grid+1)+(3,),order='F'),node_coord0(grid,size,origin),atol=atol):
raise ValueError('Input data is not ordered (x fast, z slow).')
return (grid,size,origin)
@ -376,15 +404,15 @@ def regrid(size,F,new_grid):
Parameters
----------
size : numpy.ndarray
size : numpy.ndarray of shape (3)
physical size
F : numpy.ndarray
F : numpy.ndarray of shape (:,:,:,3,3)
deformation gradient field
new_grid : numpy.ndarray
new_grid : numpy.ndarray of shape (3)
new grid for undeformed coordinates
"""
c = cell_coord0(F.shape[:3][::-1],size) \
c = cell_coord0(F.shape[:3],size) \
+ cell_displacement_avg(size,F) \
+ cell_displacement_fluct(size,F)

47
python/tests/test_grid_filters.py
View File

@ -9,13 +9,13 @@ class TestGridFilters:
size = np.random.random(3)
grid = np.random.randint(8,32,(3))
coord = grid_filters.cell_coord0(grid,size)
assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid[::-1]) + (3,)
assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid) + (3,)
def test_node_coord0(self):
size = np.random.random(3)
grid = np.random.randint(8,32,(3))
coord = grid_filters.node_coord0(grid,size)
assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid[::-1]+1) + (3,)
assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid+1) + (3,)
def test_coord0(self):
size = np.random.random(3)
@ -31,7 +31,7 @@ class TestGridFilters:
size = np.random.random(3)
origin = np.random.random(3)
coord0 = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode)) # noqa
_grid,_size,_origin = eval('grid_filters.{}_coord0_gridSizeOrigin(coord0.reshape(-1,3))'.format(mode))
_grid,_size,_origin = eval('grid_filters.{}_coord0_gridSizeOrigin(coord0.reshape(-1,3,order="F"))'.format(mode))
assert np.allclose(grid,_grid) and np.allclose(size,_size) and np.allclose(origin,_origin)
def test_displacement_fluct_equivalence(self):
@ -57,9 +57,9 @@ class TestGridFilters:
shifted = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode))
unshifted = eval('grid_filters.{}_coord0(grid,size)'.format(mode))
if mode == 'cell':
assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]) +(3,)))
assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid) +(3,)))
elif mode == 'node':
assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]+1)+(3,)))
assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid+1)+(3,)))
@pytest.mark.parametrize('function',[grid_filters.cell_displacement_avg,
grid_filters.node_displacement_avg])
@ -80,8 +80,43 @@ class TestGridFilters:
F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3))
assert np.allclose(function(size,F),0.0)
@pytest.mark.parametrize('function',[grid_filters.coord0_check,
grid_filters.node_coord0_gridSizeOrigin,
grid_filters.cell_coord0_gridSizeOrigin])
def test_invalid_coordinates(self,function):
invalid_coordinates = np.random.random((np.random.randint(12,52),3))
with pytest.raises(ValueError):
function(invalid_coordinates)
@pytest.mark.parametrize('function',[grid_filters.node_coord0_gridSizeOrigin,
grid_filters.cell_coord0_gridSizeOrigin])
def test_uneven_spaced_coordinates(self,function):
start = np.random.random(3)
end = np.random.random(3)*10. + start
grid = np.random.randint(8,32,(3))
uneven = np.stack(np.meshgrid(np.logspace(start[0],end[0],grid[0]),
np.logspace(start[1],end[1],grid[1]),
np.logspace(start[2],end[2],grid[2]),indexing = 'ij'),
axis = -1).reshape((grid.prod(),3),order='F')
with pytest.raises(ValueError):
function(uneven)
@pytest.mark.parametrize('mode',[True,False])
@pytest.mark.parametrize('function',[grid_filters.node_coord0_gridSizeOrigin,
grid_filters.cell_coord0_gridSizeOrigin])
def test_unordered_coordinates(self,function,mode):
origin = np.random.random(3)
size = np.random.random(3)*10.+origin
grid = np.random.randint(8,32,(3))
unordered = grid_filters.node_coord0(grid,size,origin).reshape(-1,3)
if mode:
with pytest.raises(ValueError):
function(unordered,mode)
else:
function(unordered,mode)
def test_regrid(self):
size = np.random.random(3)
grid = np.random.randint(8,32,(3))
F = np.broadcast_to(np.eye(3), tuple(grid[::-1])+(3,3))
F = np.broadcast_to(np.eye(3), tuple(grid)+(3,3))
assert all(grid_filters.regrid(size,F,grid) == np.arange(grid.prod()))

74
src/constitutive.f90
View File

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

2
src/constitutive_plastic_disloUCLA.f90
View File

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

2
src/constitutive_plastic_dislotwin.f90
View File

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

2
src/constitutive_plastic_isotropic.f90
View File

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

2
src/constitutive_plastic_kinehardening.f90
View File

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

2
src/constitutive_plastic_none.f90
View File

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

5
src/constitutive_plastic_nonlocal.f90
View File

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

2
src/constitutive_plastic_phenopowerlaw.f90
View File

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

544
src/crystallite.f90
View File

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

66
src/material.f90
View File

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

23
src/numerics.f90
View File

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

3
src/prec.f90
View File

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

2
src/source_damage_anisoBrittle.f90
View File

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

2
src/source_damage_anisoDuctile.f90
View File

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

2
src/source_damage_isoBrittle.f90
View File

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

2
src/source_damage_isoDuctile.f90
View File

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

2
src/source_thermal_dissipation.f90
View File

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

2
src/source_thermal_externalheat.f90
View File

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