diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index dd20c8baa..501b9e193 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -203,7 +203,6 @@ Post_OrientationConversion: stage: postprocessing script: - OrientationConversion/test.py - - OrientationConversion/test2.py except: - master - release diff --git a/VERSION b/VERSION index 7be7a7f6c..3cd762d43 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2303-g2a6132b7 +v2.0.3-2390-g524706d3 diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index 70c883aeb..2230ba563 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -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 diff --git a/processing/post/addCompatibilityMismatch.py b/processing/post/addCompatibilityMismatch.py index 0678e78ba..59eafc077 100755 --- a/processing/post/addCompatibilityMismatch.py +++ b/processing/post/addCompatibilityMismatch.py @@ -16,8 +16,8 @@ scriptID = ' '.join([scriptName,damask.version]) 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. @@ -62,19 +62,18 @@ def volTetrahedron(coords): def volumeMismatch(size,F,nodes): """ Calculates the volume mismatch. - - volume mismatch is defined as the difference between volume of reconstructed + + volume mismatch is defined as the difference between volume of reconstructed (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,47 +90,45 @@ 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): """ Routine to calculate the shape mismatch. - + shape mismatch is defined as difference between the vectors from the central point to 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 @@ -174,24 +171,24 @@ if filenames == []: filenames = [None] for name in filenames: damask.util.report(scriptName,name) - + 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) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 87d1ab2f6..17459a2df 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -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) diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py index f74d876bc..6f4a60192 100755 --- a/processing/post/addDisplacement.py +++ b/processing/post/addDisplacement.py @@ -5,8 +5,6 @@ import sys from io import StringIO from optparse import OptionParser -import numpy as np - import damask @@ -51,23 +49,23 @@ 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) diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 2619bc499..50048b44e 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -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) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 409b2ce6d..f63f24789 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -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) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index 137b8121e..e1752b0a4 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -172,7 +172,7 @@ for name in filenames: elif inputtype == 'matrix': d = representations['matrix'][1] - o = damask.Rotation.fromMatrix(list(map(float,table.data[column:column+d]))) + o = damask.Rotation.fromMatrix(np.array(list(map(float,table.data[column:column+d]))).reshape(3,3)) elif inputtype == 'frame': M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ diff --git a/processing/post/addSchmidfactors.py b/processing/post/addSchmidfactors.py index 94c062766..b3c4c7500 100755 --- a/processing/post/addSchmidfactors.py +++ b/processing/post/addSchmidfactors.py @@ -214,7 +214,7 @@ for name in filenames: outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - o = damask.Rotation(list(map(float,table.data[column:column+4]))) + o = damask.Rotation(np.array(list(map(float,table.data[column:column+4])))) table.data_append( np.abs( np.sum(slip_direction * (o * force) ,axis=1) \ * np.sum(slip_normal * (o * normal),axis=1))) diff --git a/processing/post/averageDown.py b/processing/post/averageDown.py index 39d925d19..341cc748d 100755 --- a/processing/post/averageDown.py +++ b/processing/post/averageDown.py @@ -61,7 +61,7 @@ if any(shift != 0): prefix += 'shift{:+}{:+}{:+}_'.format(*shift) for name in filenames: damask.util.report(scriptName,name) - + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) if (options.grid is None or options.size is None): @@ -87,11 +87,11 @@ for name in filenames: origin = list(-(packing//2)) + [0])\ [::packing[0],::packing[1],::packing[2],:].reshape((packedGrid.prod(),-1),order = 'F') - + 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)) diff --git a/processing/post/blowUp.py b/processing/post/blowUp.py index 316b74753..23c6b2ef2 100755 --- a/processing/post/blowUp.py +++ b/processing/post/blowUp.py @@ -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)) diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py index df40176d8..e6755c28e 100755 --- a/processing/pre/geom_fromVoronoiTessellation.py +++ b/processing/pre/geom_fromVoronoiTessellation.py @@ -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) diff --git a/processing/pre/seeds_fromGeom.py b/processing/pre/seeds_fromGeom.py index 80dc0d6f5..0b741a077 100755 --- a/processing/pre/seeds_fromGeom.py +++ b/processing/pre/seeds_fromGeom.py @@ -45,7 +45,7 @@ options.blacklist = [int(i) for i in options.blacklist] for name in filenames: damask.util.report(scriptName,name) - + geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) microstructure = geom.get_microstructure().reshape((-1,1),order='F') @@ -53,9 +53,9 @@ for name in filenames: np.full(geom.grid.prod(),True,dtype=bool), 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:]), 'grid\ta {}\tb {}\tc {}'.format(*geom.grid), diff --git a/processing/pre/seeds_fromRandom.py b/processing/pre/seeds_fromRandom.py index 2de513c2c..be690713c 100755 --- a/processing/pre/seeds_fromRandom.py +++ b/processing/pre/seeds_fromRandom.py @@ -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: diff --git a/python/.coveragerc b/python/.coveragerc index c712d2595..5daa25bb2 100644 --- a/python/.coveragerc +++ b/python/.coveragerc @@ -1,2 +1,5 @@ [run] omit = tests/* + damask/_asciitable.py + damask/_test.py + damask/config/* diff --git a/python/damask/_geom.py b/python/damask/_geom.py index e2c2428fe..892000b7c 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -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) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 76475e057..e6c732007 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -38,6 +38,9 @@ class Orientation: else: self.rotation = Rotation.fromQuaternion(rotation) # assume quaternion + if self.rotation.quaternion.shape != (4,): + raise NotImplementedError('Support for multiple rotations missing') + def disorientation(self, other, SST = True, diff --git a/python/damask/_result.py b/python/damask/_result.py index a8181b7a1..9171682ff 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -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'][()] diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index a2ea47007..716977a79 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -1,6 +1,7 @@ import numpy as np from ._Lambert import ball_to_cube, cube_to_ball +from . import mechanics _P = -1 @@ -61,6 +62,8 @@ class Rotation: def __repr__(self): """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles.""" + if self.quaternion.shape != (4,): + raise NotImplementedError('Support for multiple rotations missing') return '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), 'Matrix:\n{}'.format(self.asMatrix()), @@ -83,6 +86,8 @@ class Rotation: considere rotation of (3,3,3,3)-matrix """ + if self.quaternion.shape != (4,): + raise NotImplementedError('Support for multiple rotations missing') if isinstance(other, Rotation): # rotate a rotation self_q = self.quaternion[0] self_p = self.quaternion[1:] @@ -107,7 +112,7 @@ class Rotation: elif other.shape == (3,3,): # rotate a single (3x3)-matrix return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) elif other.shape == (3,3,3,3,): - raise NotImplementedError + raise NotImplementedError('Support for rotation of 4th order tensors missing') else: return NotImplemented else: @@ -116,7 +121,7 @@ class Rotation: def inverse(self): """In-place inverse rotation/backward rotation.""" - self.quaternion[1:] *= -1 + self.quaternion[...,1:] *= -1 return self def inversed(self): @@ -125,12 +130,12 @@ class Rotation: def standardize(self): - """In-place quaternion representation with positive q.""" - if self.quaternion[0] < 0.0: self.quaternion*=-1 + """In-place quaternion representation with positive real part.""" + self.quaternion[self.quaternion[...,0] < 0.0] *= -1 return self def standardized(self): - """Quaternion representation with positive q.""" + """Quaternion representation with positive real part.""" return self.copy().standardize() @@ -157,15 +162,17 @@ class Rotation: Rotation from which the average is rotated. """ + if self.quaternion.shape != (4,) or other.quaternion.shape != (4,): + raise NotImplementedError('Support for multiple rotations missing') return Rotation.fromAverage([self,other]) ################################################################################################ # convert to different orientation representations (numpy arrays) - def asQuaternion(self): + def as_quaternion(self): """ - Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object. + Unit quaternion [q, p_1, p_2, p_3]. Parameters ---------- @@ -175,8 +182,8 @@ class Rotation: """ return self.quaternion - def asEulers(self, - degrees = False): + def as_Eulers(self, + degrees = False): """ Bunge-Euler angles: (φ_1, ϕ, φ_2). @@ -190,9 +197,9 @@ class Rotation: if degrees: eu = np.degrees(eu) return eu - def asAxisAngle(self, - degrees = False, - pair = False): + def as_axis_angle(self, + degrees = False, + pair = False): """ Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω). @@ -205,15 +212,15 @@ class Rotation: """ ax = Rotation.qu2ax(self.quaternion) - if degrees: ax[3] = np.degrees(ax[3]) - return (ax[:3],ax[3]) if pair else ax + if degrees: ax[...,3] = np.degrees(ax[...,3]) + return (ax[...,:3],ax[...,3]) if pair else ax - def asMatrix(self): + def as_matrix(self): """Rotation matrix.""" return Rotation.qu2om(self.quaternion) - def asRodrigues(self, - vector = False): + def as_Rodrigues(self, + vector = False): """ Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2). @@ -224,9 +231,9 @@ class Rotation: """ ro = Rotation.qu2ro(self.quaternion) - return ro[:3]*ro[3] if vector else ro + return ro[...,:3]*ro[...,3] if vector else ro - def asHomochoric(self): + def as_homochoric(self): """Homochoric vector: (h_1, h_2, h_3).""" return Rotation.qu2ho(self.quaternion) @@ -234,7 +241,7 @@ class Rotation: """Cubochoric vector: (c_1, c_2, c_3).""" return Rotation.qu2cu(self.quaternion) - def asM(self): + def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M """ Intermediate representation supporting quaternion averaging. @@ -244,114 +251,133 @@ class Rotation: https://doi.org/10.2514/1.28949 """ - return np.outer(self.quaternion,self.quaternion) + return np.einsum('...i,...j',self.quaternion,self.quaternion) + # for compatibility (old names do not follow convention) + asM = M + asQuaternion = as_quaternion + asEulers = as_Eulers + asAxisAngle = as_axis_angle + asMatrix = as_matrix + asRodrigues = as_Rodrigues + asHomochoric = as_homochoric ################################################################################################ - # static constructors. The input data needs to follow the convention, options allow to - # relax these convections + # Static constructors. The input data needs to follow the conventions, options allow to + # relax the conventions. @staticmethod - def fromQuaternion(quaternion, - acceptHomomorph = False, - P = -1): + def from_quaternion(quaternion, + acceptHomomorph = False, + P = -1): - qu = quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \ - else np.array(quaternion,dtype=float) - if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1 - if qu[0] < 0.0: - if acceptHomomorph: - qu *= -1. - else: - raise ValueError('Quaternion has negative first component: {}.'.format(qu[0])) - if not np.isclose(np.linalg.norm(qu), 1.0): - raise ValueError('Quaternion is not of unit length: {} {} {} {}.'.format(*qu)) + qu = np.array(quaternion,dtype=float) + if qu.shape[:-2:-1] != (4,): + raise ValueError('Invalid shape.') + + if P > 0: qu[...,1:4] *= -1 # convert from P=1 to P=-1 + if acceptHomomorph: + qu[qu[...,0] < 0.0] *= -1 + else: + if np.any(qu[...,0] < 0.0): + raise ValueError('Quaternion with negative first (real) component.') + if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0)): + raise ValueError('Quaternion is not of unit length.') return Rotation(qu) @staticmethod - def fromEulers(eulers, - degrees = False): + def from_Eulers(eulers, + degrees = False): + + eu = np.array(eulers,dtype=float) + if eu.shape[:-2:-1] != (3,): + raise ValueError('Invalid shape.') - eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \ - else np.array(eulers,dtype=float) eu = np.radians(eu) if degrees else eu - if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi: - raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]: {} {} {}.'.format(*eu)) + if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or np.any(eu[...,1] > np.pi): # ToDo: No separate check for PHI + raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].') return Rotation(Rotation.eu2qu(eu)) @staticmethod - def fromAxisAngle(angleAxis, - degrees = False, - normalise = False, - P = -1): + def from_axis_angle(axis_angle, + degrees = False, + normalise = False, + P = -1): - ax = angleAxis if isinstance(angleAxis, np.ndarray) and angleAxis.dtype == np.dtype(float) \ - else np.array(angleAxis,dtype=float) - if P > 0: ax[0:3] *= -1 # convert from P=1 to P=-1 - if degrees: ax[ 3] = np.radians(ax[3]) - if normalise: ax[0:3] /= np.linalg.norm(ax[0:3]) - if ax[3] < 0.0 or ax[3] > np.pi: - raise ValueError('Axis angle rotation angle outside of [0..π]: {}.'.format(ax[3])) - if not np.isclose(np.linalg.norm(ax[0:3]), 1.0): - raise ValueError('Axis angle rotation axis is not of unit length: {} {} {}.'.format(*ax[0:3])) + ax = np.array(axis_angle,dtype=float) + if ax.shape[:-2:-1] != (4,): + raise ValueError('Invalid shape.') + + if P > 0: ax[...,0:3] *= -1 # convert from P=1 to P=-1 + if degrees: ax[..., 3] = np.radians(ax[...,3]) + if normalise: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1) + if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi): + raise ValueError('Axis angle rotation angle outside of [0..π].') + if not np.all(np.isclose(np.linalg.norm(ax[...,0:3],axis=-1), 1.0)): + raise ValueError('Axis angle rotation axis is not of unit length.') return Rotation(Rotation.ax2qu(ax)) @staticmethod - def fromBasis(basis, - orthonormal = True, - reciprocal = False, - ): + def from_basis(basis, + orthonormal = True, + reciprocal = False): + + om = np.array(basis,dtype=float) + if om.shape[:-3:-1] != (3,3): + raise ValueError('Invalid shape.') - om = basis if isinstance(basis, np.ndarray) else np.array(basis).reshape(3,3) if reciprocal: - om = np.linalg.inv(om.T/np.pi) # transform reciprocal basis set + om = np.linalg.inv(mechanics.transpose(om)/np.pi) # transform reciprocal basis set orthonormal = False # contains stretch if not orthonormal: (U,S,Vh) = np.linalg.svd(om) # singular value decomposition - om = np.dot(U,Vh) - if not np.isclose(np.linalg.det(om),1.0): - raise ValueError('matrix is not a proper rotation: {}.'.format(om)) - if not np.isclose(np.dot(om[0],om[1]), 0.0) \ - or not np.isclose(np.dot(om[1],om[2]), 0.0) \ - or not np.isclose(np.dot(om[2],om[0]), 0.0): - raise ValueError('matrix is not orthogonal: {}.'.format(om)) + om = np.einsum('...ij,...jl->...il',U,Vh) + if not np.all(np.isclose(np.linalg.det(om),1.0)): + raise ValueError('Orientation matrix has determinant ≠ 1.') + if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \ + or not np.all(np.isclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0)) \ + or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)): + raise ValueError('Orientation matrix is not orthogonal.') return Rotation(Rotation.om2qu(om)) @staticmethod - def fromMatrix(om, - ): + def from_matrix(om): - return Rotation.fromBasis(om) + return Rotation.from_basis(om) @staticmethod - def fromRodrigues(rodrigues, - normalise = False, - P = -1): + def from_Rodrigues(rodrigues, + normalise = False, + P = -1): - ro = rodrigues if isinstance(rodrigues, np.ndarray) and rodrigues.dtype == np.dtype(float) \ - else np.array(rodrigues,dtype=float) - if P > 0: ro[0:3] *= -1 # convert from P=1 to P=-1 - if normalise: ro[0:3] /= np.linalg.norm(ro[0:3]) - if not np.isclose(np.linalg.norm(ro[0:3]), 1.0): - raise ValueError('Rodrigues rotation axis is not of unit length: {} {} {}.'.format(*ro[0:3])) - if ro[3] < 0.0: - raise ValueError('Rodrigues rotation angle not positive: {}.'.format(ro[3])) + ro = np.array(rodrigues,dtype=float) + if ro.shape[:-2:-1] != (4,): + raise ValueError('Invalid shape.') + + if P > 0: ro[...,0:3] *= -1 # convert from P=1 to P=-1 + if normalise: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1) + if np.any(ro[...,3] < 0.0): + raise ValueError('Rodrigues vector rotation angle not positive.') + if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)): + raise ValueError('Rodrigues vector rotation axis is not of unit length.') return Rotation(Rotation.ro2qu(ro)) @staticmethod - def fromHomochoric(homochoric, - P = -1): + def from_homochoric(homochoric, + P = -1): + + ho = np.array(homochoric,dtype=float) + if ho.shape[:-2:-1] != (3,): + raise ValueError('Invalid shape.') - ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \ - else np.array(homochoric,dtype=float) if P > 0: ho *= -1 # convert from P=1 to P=-1 - if np.linalg.norm(ho) > (3.*np.pi/4.)**(1./3.)+1e-9: - raise ValueError('Coordinate outside of the sphere: {} {} {}.'.format(ho)) + if np.any(np.linalg.norm(ho,axis=-1) > (3.*np.pi/4.)**(1./3.)+1e-9): + raise ValueError('Homochoric coordinate outside of the sphere.') return Rotation(Rotation.ho2qu(ho)) @@ -359,11 +385,12 @@ class Rotation: def fromCubochoric(cubochoric, P = -1): - cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \ - else np.array(cubochoric,dtype=float) + cu = np.array(cubochoric,dtype=float) + if cu.shape[:-2:-1] != (3,): + raise ValueError('Invalid shape.') if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9: - raise ValueError('Coordinate outside of the cube: {} {} {}.'.format(*cu)) + raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu)) ho = Rotation.cu2ho(cu) if P > 0: ho *= -1 # convert from P=1 to P=-1 @@ -403,17 +430,34 @@ class Rotation: return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) - @staticmethod - def fromRandom(): - r = np.random.random(3) - A = np.sqrt(r[2]) - B = np.sqrt(1.0-r[2]) - return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A, - np.sin(2.0*np.pi*r[1])*B, - np.cos(2.0*np.pi*r[1])*B, - np.sin(2.0*np.pi*r[0])*A])).standardize() + def from_random(shape=None): + if shape is None: + r = np.random.random(3) + elif hasattr(shape, '__iter__'): + r = np.random.random(tuple(shape)+(3,)) + else: + r = np.random.random((shape,3)) + A = np.sqrt(r[...,2]) + B = np.sqrt(1.0-r[...,2]) + q = np.stack([np.cos(2.0*np.pi*r[...,0])*A, + np.sin(2.0*np.pi*r[...,1])*B, + np.cos(2.0*np.pi*r[...,1])*B, + np.sin(2.0*np.pi*r[...,0])*A],axis=-1) + + return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q).standardize() + + + # for compatibility (old names do not follow convention) + fromQuaternion = from_quaternion + fromEulers = from_Eulers + fromAxisAngle = from_axis_angle + fromBasis = from_basis + fromMatrix = from_matrix + fromRodrigues = from_Rodrigues + fromHomochoric = from_homochoric + fromRandom = from_random #################################################################################################### # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations @@ -808,12 +852,11 @@ class Rotation: c = np.cos(ax[3]*0.5) s = np.sin(ax[3]*0.5) qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) - return qu else: c = np.cos(ax[...,3:4]*.5) s = np.sin(ax[...,3:4]*.5) qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s])) - return qu + return qu @staticmethod def ax2om(ax): @@ -859,7 +902,7 @@ class Rotation: # 180 degree case ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ [np.tan(ax[3]*0.5)] - return np.array(ro) + ro = np.array(ro) else: ro = np.block([ax[...,:3], np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0), @@ -867,7 +910,7 @@ class Rotation: np.tan(ax[...,3:4]*0.5)) ]) ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0] - return ro + return ro @staticmethod def ax2ho(ax): @@ -875,11 +918,10 @@ class Rotation: if len(ax.shape) == 1: f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) ho = ax[0:3] * f - return ho else: f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0) ho = ax[...,:3] * f - return ho + return ho @staticmethod def ax2cu(ax): @@ -936,7 +978,6 @@ class Rotation: f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi) ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape), np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0)) - return ho @staticmethod @@ -1010,7 +1051,7 @@ class Rotation: if len(ho.shape) == 1: return ball_to_cube(ho) else: - raise NotImplementedError + raise NotImplementedError('Support for multiple rotations missing') #---------- Cubochoric ---------- @@ -1045,4 +1086,4 @@ class Rotation: if len(cu.shape) == 1: return cube_to_ball(cu) else: - raise NotImplementedError + raise NotImplementedError('Support for multiple rotations missing') diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index d8b136a6b..eaee42924 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -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,13 +238,13 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True): start = origin + delta*.5 end = origin - delta*.5 + size - if not _np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0])) and \ - _np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1])) and \ - _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2])): + if not (_np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0])) and \ + _np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1])) and \ + _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2]))): raise ValueError('Regular grid spacing violated.') - if ordered and not _np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin)): - 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)): + raise ValueError('Input data is not ordered (x fast, z slow).') return (grid,size,origin) @@ -241,17 +268,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): @@ -260,7 +288,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. @@ -275,14 +303,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): @@ -291,7 +319,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. @@ -306,15 +334,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): @@ -335,14 +363,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). @@ -357,13 +385,13 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False): if (grid+1).prod() != len(coord0): raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) - if not _np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \ - _np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \ - _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1)): + if not (_np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \ + _np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \ + _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1))): raise ValueError('Regular grid spacing violated.') - if ordered and not _np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)): - 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)): + raise ValueError('Input data is not ordered (x fast, z slow).') return (grid,size,origin) @@ -374,15 +402,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) diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index e19f140fb..c81399d94 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -135,16 +135,16 @@ def PK2(P,F): Parameters ---------- - P : numpy.ndarray of shape (:,3,3) or (3,3) + P : numpy.ndarray of shape (...,3,3) or (3,3) First Piola-Kirchhoff stress. - F : numpy.ndarray of shape (:,3,3) or (3,3) + F : numpy.ndarray of shape (...,3,3) or (3,3) Deformation gradient. """ if _np.shape(F) == _np.shape(P) == (3,3): S = _np.dot(_np.linalg.inv(F),P) else: - S = _np.einsum('ijk,ikl->ijl',_np.linalg.inv(F),P) + S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P) return symmetric(S) @@ -241,7 +241,7 @@ def symmetric(T): Parameters ---------- - T : numpy.ndarray of shape (:,3,3) or (3,3) + T : numpy.ndarray of shape (...,3,3) or (3,3) Tensor of which the symmetrized values are computed. """ @@ -254,12 +254,12 @@ def transpose(T): Parameters ---------- - T : numpy.ndarray of shape (:,3,3) or (3,3) + T : numpy.ndarray of shape (...,3,3) or (3,3) Tensor of which the transpose is computed. """ return T.T if _np.shape(T) == (3,3) else \ - _np.transpose(T,(0,2,1)) + _np.swapaxes(T,axis2=-2,axis1=-1) def _polar_decomposition(T,requested): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 2ac819f4c..302e895ee 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -157,6 +157,30 @@ class TestRotation: print(m,o,rot.asQuaternion()) assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 + @pytest.mark.parametrize('function',[Rotation.from_quaternion, + Rotation.from_Eulers, + Rotation.from_axis_angle, + Rotation.from_matrix, + Rotation.from_Rodrigues, + Rotation.from_homochoric]) + def test_invalid_shape(self,function): + invalid_shape = np.random.random(np.random.randint(8,32,(3))) + with pytest.raises(ValueError): + function(invalid_shape) + + @pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])), + (Rotation.from_quaternion, np.array([1,1,1,0])), + (Rotation.from_Eulers, np.array([1,4,0])), + (Rotation.from_axis_angle, np.array([1,0,0,4])), + (Rotation.from_axis_angle, np.array([1,1,0,1])), + (Rotation.from_matrix, np.random.rand(3,3)), + (Rotation.from_Rodrigues, np.array([1,0,0,-1])), + (Rotation.from_Rodrigues, np.array([1,1,0,1])), + (Rotation.from_homochoric, np.array([2,2,2])) ]) + def test_invalid(self,function,invalid): + with pytest.raises(ValueError): + function(invalid) + @pytest.mark.parametrize('conversion',[Rotation.qu2om, Rotation.qu2eu, Rotation.qu2ax, diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index acbdbf688..eb359006a 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -4,18 +4,18 @@ import numpy as np from damask import grid_filters class TestGridFilters: - + def test_cell_coord0(self): 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())) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index c590a86b5..5394fc90c 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -73,28 +73,24 @@ subroutine CPFEM_initAll(el,ip) integer(pInt), intent(in) :: el, & !< FE el number ip !< FE integration point number - !$OMP CRITICAL(init) - if (.not. CPFEM_init_done) then - call DAMASK_interface_init - call prec_init - call IO_init - call numerics_init - call debug_init - call config_init - call math_init - call rotations_init - call HDF5_utilities_init - call results_init - call discretization_marc_init(ip, el) - call lattice_init - call material_init - call constitutive_init - call crystallite_init - call homogenization_init - call CPFEM_init - CPFEM_init_done = .true. - endif - !$OMP END CRITICAL(init) + CPFEM_init_done = .true. + call DAMASK_interface_init + call prec_init + call IO_init + call numerics_init + call debug_init + call config_init + call math_init + call rotations_init + call HDF5_utilities_init + call results_init + call discretization_marc_init(ip, el) + call lattice_init + call material_init + call constitutive_init + call crystallite_init + call homogenization_init + call CPFEM_init end subroutine CPFEM_initAll diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 9fd41459b..fa0711b02 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -261,11 +261,10 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc + !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS if (.not. CPFEM_init_done) call CPFEM_initAll(m(1),nn) - !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS - computationMode = 0 ! save initialization value, since it does not result in any calculation if (lovl == 4 ) then ! jacobian requested by marc if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 94e590ab3..e2c9dbc05 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -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 !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 index 6be86f266..90a933910 100644 --- a/src/constitutive_plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -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 diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 21fe555f2..7c7d24ab8 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -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 diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index 94fc9817d..ecf029124 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -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 diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index 5843f5b5e..36b1eedf9 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -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 diff --git a/src/constitutive_plastic_none.f90 b/src/constitutive_plastic_none.f90 index 7ff1c76f7..667fe5638 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -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 diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index f872e3846..77f1556f6 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -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,:) diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index a980d6106..12a30478a 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -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 diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 9bc254e0c..929fec862 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -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,7 +301,9 @@ 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 + todo = .false. #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & .and. FEsolving_execElem(1) <= debug_e & @@ -344,7 +346,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 +363,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 +382,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 +417,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 +440,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 +612,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 +781,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 +838,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 +851,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 +984,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 +991,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 +1000,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 +1011,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 +1146,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 +1159,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 +1190,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 +1211,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 +1224,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), & - 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 + 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) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1274,36 +1261,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), & - 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 + 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) 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 +1285,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 +1293,7 @@ subroutine integrateStateAdaptiveEuler enddo endif + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1328,189 +1302,70 @@ subroutine integrateStateAdaptiveEuler end subroutine integrateStateAdaptiveEuler -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with 4th order explicit Runge Kutta method -!-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4 +!--------------------------------------------------------------------------------------------------- +!> @brief Integrate state (including stress integration) with the classic Runge Kutta method +!--------------------------------------------------------------------------------------------------- +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) + 0.0_pReal, 0.0_pReal, 1.0_pReal],& + shape(A)) + real(pReal), dimension(3), parameter :: & + C = [0.5_pReal, 0.5_pReal, 1.0_pReal] + 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] - integer :: & - e, & ! element index in element loop - i, & ! integration point index in ip loop - g, & ! grain index in grain loop - stage, & ! stage index in integration stage loop - n, & - p, & - c, & - s, & - sizeDotState - logical :: & - nonlocalBroken - - real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState - real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState - nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RK4dotState,source_RK4dotState) - 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); c = material_phaseMemberAt(g,i,e) - - 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) .or. crystallite_localPlasticity(g,i,e))) & - nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle - - do stage = 1,3 - sizeDotState = plasticState(p)%sizeDotState - plastic_RK4dotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RK4dotState(1:sizeDotState,1) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - source_RK4dotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RK4dotState(1:sizeDotState,1,s) - enddo - - do n = 2, stage - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + A(n,stage) * plastic_RK4dotState(1:sizeDotState,n) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & - + A(n,stage) * source_RK4dotState(1:sizeDotState,n,s) - enddo - enddo - - 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) - 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) - 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) - - 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), & - 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 - - enddo - - if(.not. crystallite_todo(g,i,e)) cycle - - sizeDotState = plasticState(p)%sizeDotState - - plastic_RK4dotState(1:sizeDotState,4) = plasticState (p)%dotState(:,c) - - plasticState(p)%dotState(:,c) = matmul(plastic_RK4dotState(1:sizeDotState,1:4),B) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - - source_RK4dotState(1:sizeDotState,4,s) = sourceState(p)%p(s)%dotState(:,c) - - sourceState(p)%p(s)%dotState(:,c) = matmul(source_RK4dotState(1:sizeDotState,1:4,s),B) - 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) - 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) - - 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 - - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - - if (nonlocalBroken) call nonlocalConvergenceCheck + call integrateStateRK(todo,A,B,C) 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 +!--------------------------------------------------------------------------------------------------- +!> @brief Integrate state (including stress integration) with the Cash-Carp method +!--------------------------------------------------------------------------------------------------- +subroutine integrateStateRKCK45(todo) + + logical, dimension(:,:,:), intent(in) :: todo real(pReal), dimension(5,5), parameter :: & A = reshape([& - .2_pReal, .075_pReal, .3_pReal, -11.0_pReal/54.0_pReal, 1631.0_pReal/55296.0_pReal, & - .0_pReal, .225_pReal, -.9_pReal, 2.5_pReal, 175.0_pReal/512.0_pReal, & - .0_pReal, .0_pReal, 1.2_pReal, -70.0_pReal/27.0_pReal, 575.0_pReal/13824.0_pReal, & - .0_pReal, .0_pReal, .0_pReal, 35.0_pReal/27.0_pReal, 44275.0_pReal/110592.0_pReal, & - .0_pReal, .0_pReal, .0_pReal, .0_pReal, 253.0_pReal/4096.0_pReal], & - [5,5], order=[2,1]) !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6) - + 1._pReal/5._pReal, .0_pReal, .0_pReal, .0_pReal, .0_pReal, & + 3._pReal/40._pReal, 9._pReal/40._pReal, .0_pReal, .0_pReal, .0_pReal, & + 3_pReal/10._pReal, -9._pReal/10._pReal, 6._pReal/5._pReal, .0_pReal, .0_pReal, & + -11._pReal/54._pReal, 5._pReal/2._pReal, -70.0_pReal/27.0_pReal, 35.0_pReal/27.0_pReal, .0_pReal, & + 1631._pReal/55296._pReal,175._pReal/512._pReal,575._pReal/13824._pReal,44275._pReal/110592._pReal,253._pReal/4096._pReal],& + shape(A)) + real(pReal), dimension(5), parameter :: & + C = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] real(pReal), dimension(6), parameter :: & B = & [37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, & - 125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], & !< coefficients in Butcher tableau (used for final integration and error estimate) + 125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], & DB = B - & [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& - 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 0.25_pReal] !< coefficients in Butcher tableau (used for final integration and error estimate) + 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - real(pReal), dimension(5), parameter :: & - CC = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6) + call integrateStateRK(todo,A,B,C,DB) + +end subroutine integrateStateRKCK45 + + +!-------------------------------------------------------------------------------------------------- +!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an +!! embedded explicit Runge-Kutta method +!-------------------------------------------------------------------------------------------------- +subroutine integrateStateRK(todo,A,B,CC,DB) + + logical, dimension(:,:,:), intent(in) :: todo + + real(pReal), dimension(:,:), intent(in) :: A + real(pReal), dimension(:), intent(in) :: B, CC + real(pReal), dimension(:), intent(in), optional :: DB integer :: & e, & ! element index in element loop @@ -1523,33 +1378,29 @@ subroutine integrateStateRKCK45 s, & sizeDotState logical :: & - nonlocalBroken - real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState - real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState + nonlocalBroken, broken + real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: 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), & - 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 + 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) cycle - do stage = 1,5 + do stage = 1,size(A,1) sizeDotState = plasticState(p)%sizeDotState plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) @@ -1581,83 +1432,64 @@ 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), & - 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 + 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,p,c) + if(broken) exit enddo - - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState - plastic_RKdotState(1:sizeDotState,6) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:6),B) + plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) + plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) 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) & - * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + if(present(DB)) & + broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & + * crystallite_subdt(g,i,e), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - source_RKdotState(1:sizeDotState,6,s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:6,s),B) + source_RKdotState(1:sizeDotState,size(B),s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:size(B),s),B) 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. & - 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)) + if(present(DB)) & + broken = broken .or. .not. converged(matmul(source_RKdotState(1:sizeDotState,1:size(DB),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 +end subroutine integrateStateRK !-------------------------------------------------------------------------------------------------- @@ -1666,7 +1498,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 +1529,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 diff --git a/src/material.f90 b/src/material.f90 index aefe44878..749c9a3d8 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -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 diff --git a/src/numerics.f90 b/src/numerics.f90 index a29601322..8d242c71d 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -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') diff --git a/src/prec.f90 b/src/prec.f90 index dd39ddc05..646f7dd69 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -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 diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 3978be959..b3af24f38 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -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' diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 2b818e2cf..79cc0c2f7 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -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' diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index ed6d89a89..9eacb4516 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -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' diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 7024e595a..96754725d 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -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' diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 521c79077..c323e68b5 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -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 diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index ade13bef2..06b8a5197 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -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