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

Merge remote-tracking branch 'origin/development' into seeds-module

This commit is contained in:
Martin Diehl 2020-09-24 20:56:32 +02:00
parent 9bce1ddc29 28cded3669
commit b5ea04424b
62 changed files with 580 additions and 615 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 86f77da4aec6cb52bbcc2724f00a6cf6a7dc6e91 Subproject commit b73dcfe746eadce92ed0ab08a3bfbb19f5c8eb0a

2
VERSION
View File

@ -1 +1 @@
v3.0.0-alpha-275-g7801f527f v3.0.0-alpha-321-g0f6495430

40
processing/pre/geom_grainGrowth.py
View File

@ -66,13 +66,13 @@ for name in filenames:
grid_original = geom.grid grid_original = geom.grid
damask.util.croak(geom) damask.util.croak(geom)
materials = np.tile(geom.materials,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1 material = np.tile(geom.material,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1
grid = np.array(materials.shape) grid = np.array(material.shape)
# --- initialize support data --------------------------------------------------------------------- # --- initialize support data ---------------------------------------------------------------------
# store a copy of the initial material indices to find locations of immutable indices # store a copy of the initial material indices to find locations of immutable indices
materials_original = np.copy(materials) material_original = np.copy(material)
if not options.ndimage: if not options.ndimage:
X,Y,Z = np.mgrid[0:grid[0],0:grid[1],0:grid[2]] X,Y,Z = np.mgrid[0:grid[0],0:grid[1],0:grid[2]]
@ -88,14 +88,14 @@ for name in filenames:
for smoothIter in range(options.N): for smoothIter in range(options.N):
interfaceEnergy = np.zeros(materials.shape,dtype=np.float32) interfaceEnergy = np.zeros(material.shape,dtype=np.float32)
for i in (-1,0,1): for i in (-1,0,1):
for j in (-1,0,1): for j in (-1,0,1):
for k in (-1,0,1): for k in (-1,0,1):
# assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood) # assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood)
interfaceEnergy = np.maximum(interfaceEnergy, interfaceEnergy = np.maximum(interfaceEnergy,
getInterfaceEnergy(materials,np.roll(np.roll(np.roll( getInterfaceEnergy(material,np.roll(np.roll(np.roll(
materials,i,axis=0), j,axis=1), k,axis=2))) material,i,axis=0), j,axis=1), k,axis=2)))
# periodically extend interfacial energy array by half a grid size in positive and negative directions # periodically extend interfacial energy array by half a grid size in positive and negative directions
periodic_interfaceEnergy = np.tile(interfaceEnergy,(3,3,3))[grid[0]//2:-grid[0]//2, periodic_interfaceEnergy = np.tile(interfaceEnergy,(3,3,3))[grid[0]//2:-grid[0]//2,
@ -143,33 +143,33 @@ for name in filenames:
return_distances = False, return_distances = False,
return_indices = True) # want index of closest bulk grain return_indices = True) # want index of closest bulk grain
periodic_materials = np.tile(materials,(3,3,3))[grid[0]//2:-grid[0]//2, periodic_material = np.tile(material,(3,3,3))[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2, grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # periodically extend the geometry grid[2]//2:-grid[2]//2] # periodically extend the geometry
materials = periodic_materials[index[0], material = periodic_material[index[0],
index[1], index[1],
index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2, index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2, grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # extent grains into interface region grid[2]//2:-grid[2]//2] # extent grains into interface region
# replace immutable materials with closest mutable ones # replace immutable materials with closest mutable ones
index = ndimage.morphology.distance_transform_edt(np.in1d(materials,options.immutable).reshape(grid), index = ndimage.morphology.distance_transform_edt(np.in1d(material,options.immutable).reshape(grid),
return_distances = False, return_distances = False,
return_indices = True) return_indices = True)
materials = materials[index[0], material = material[index[0],
index[1], index[1],
index[2]] index[2]]
immutable = np.zeros(materials.shape, dtype=np.bool) immutable = np.zeros(material.shape, dtype=np.bool)
# find locations where immutable materials have been in original structure # find locations where immutable materials have been in original structure
for micro in options.immutable: for micro in options.immutable:
immutable += materials_original == micro immutable += material_original == micro
# undo any changes involving immutable materials # undo any changes involving immutable materials
materials = np.where(immutable, materials_original,materials) material = np.where(immutable, material_original,material)
damask.Geom(materials = materials[0:grid_original[0],0:grid_original[1],0:grid_original[2]], damask.Geom(material = material[0:grid_original[0],0:grid_original[1],0:grid_original[2]],
size = geom.size, size = geom.size,
origin = geom.origin, origin = geom.origin,
comments = geom.comments + [scriptID + ' ' + ' '.join(sys.argv[1:])], comments = geom.comments + [scriptID + ' ' + ' '.join(sys.argv[1:])],

12
processing/pre/mentat_spectralBox.py
View File

@ -100,7 +100,7 @@ def mesh(r,d):
#------------------------------------------------------------------------------------------------- #-------------------------------------------------------------------------------------------------
def material(): def materials():
return [\ return [\
"*new_mater standard", "*new_mater standard",
"*mater_option general:state:solid", "*mater_option general:state:solid",
@ -130,10 +130,10 @@ def geometry():
#------------------------------------------------------------------------------------------------- #-------------------------------------------------------------------------------------------------
def initial_conditions(materials): def initial_conditions(material):
elements = [] elements = []
element = 0 element = 0
for id in materials: for id in material:
element += 1 element += 1
if len(elements) < id: if len(elements) < id:
for i in range(id-len(elements)): for i in range(id-len(elements)):
@ -197,14 +197,14 @@ for name in filenames:
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
materials = geom.materials.flatten(order='F') material = geom.material.flatten(order='F')
cmds = [\ cmds = [\
init(), init(),
mesh(geom.grid,geom.size), mesh(geom.grid,geom.size),
material(), materials(),
geometry(), geometry(),
initial_conditions(materials), initial_conditions(material),
'*identify_sets', '*identify_sets',
'*show_model', '*show_model',
'*redraw', '*redraw',

14
processing/pre/seeds_fromDistribution.py
View File

@ -101,8 +101,8 @@ class myThread (threading.Thread):
#--- evaluate current seeds file ------------------------------------------------------------------ #--- evaluate current seeds file ------------------------------------------------------------------
perturbedGeom = damask.Geom.load_ASCII(perturbedGeomVFile) perturbedGeom = damask.Geom.load_ASCII(perturbedGeomVFile)
myNmaterials = len(np.unique(perturbedGeom.materials)) myNmaterials = len(np.unique(perturbedGeom.material))
currentData=np.bincount(perturbedGeom.materials.ravel())[1:]/points currentData = np.bincount(perturbedGeom.material.ravel())[1:]/points
currentError=[] currentError=[]
currentHist=[] currentHist=[]
for i in range(nMaterials): # calculate the deviation in all bins per histogram for i in range(nMaterials): # calculate the deviation in all bins per histogram
@ -217,8 +217,8 @@ points = np.array(options.grid).prod().astype('float')
# ----------- calculate target distribution and bin edges # ----------- calculate target distribution and bin edges
targetGeom = damask.Geom.load_ASCII(os.path.splitext(os.path.basename(options.target))[0]+'.geom') targetGeom = damask.Geom.load_ASCII(os.path.splitext(os.path.basename(options.target))[0]+'.geom')
nMaterials = len(np.unique(targetGeom.materials)) nMaterials = len(np.unique(targetGeom.material))
targetVolFrac = np.bincount(targetGeom.materials.flatten())/targetGeom.grid.prod().astype(np.float) targetVolFrac = np.bincount(targetGeom.material.flatten())/targetGeom.grid.prod().astype(np.float)
target = [] target = []
for i in range(1,nMaterials+1): for i in range(1,nMaterials+1):
targetHist,targetBins = np.histogram(targetVolFrac,bins=i) #bin boundaries targetHist,targetBins = np.histogram(targetVolFrac,bins=i) #bin boundaries
@ -244,10 +244,10 @@ initialGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+
initialGeomVFile.seek(0) initialGeomVFile.seek(0)
initialGeom = damask.Geom.load_ASCII(initialGeomVFile) initialGeom = damask.Geom.load_ASCII(initialGeomVFile)
if len(np.unique(targetGeom.materials)) != nMaterials: if len(np.unique(targetGeom.material)) != nMaterials:
damask.util.croak('error. Material count mismatch') damask.util.croak('error. Material count mismatch')
initialData = np.bincount(initialGeom.materials.flatten())/points initialData = np.bincount(initialGeom.material.flatten())/points
for i in range(nMaterials): for i in range(nMaterials):
initialHist = np.histogram(initialData,bins=target[i]['bins'])[0] initialHist = np.histogram(initialData,bins=target[i]['bins'])[0]
target[i]['error']=np.sqrt(np.square(np.array(target[i]['histogram']-initialHist)).sum()) target[i]['error']=np.sqrt(np.square(np.array(target[i]['histogram']-initialHist)).sum())
@ -263,7 +263,7 @@ for i in range(nMaterials):
if options.maxseeds < 1: if options.maxseeds < 1:
maxSeeds = len(np.unique(initialGeom.materials)) maxSeeds = len(np.unique(initialGeom.material))
else: else:
maxSeeds = options.maxseeds maxSeeds = options.maxseeds

8
processing/pre/seeds_fromGeom.py
View File

@ -47,11 +47,11 @@ for name in filenames:
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
materials = geom.materials.reshape((-1,1),order='F') material = geom.material.reshape((-1,1),order='F')
mask = np.logical_and(np.in1d(materials,options.whitelist,invert=False) if options.whitelist else \ mask = np.logical_and(np.in1d(material,options.whitelist,invert=False) if options.whitelist else \
np.full(geom.grid.prod(),True,dtype=bool), np.full(geom.grid.prod(),True,dtype=bool),
np.in1d(materials,options.blacklist,invert=True) if options.blacklist else \ np.in1d(material,options.blacklist,invert=True) if options.blacklist else \
np.full(geom.grid.prod(),True,dtype=bool)) np.full(geom.grid.prod(),True,dtype=bool))
seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F') seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
@ -64,5 +64,5 @@ for name in filenames:
] ]
damask.Table(seeds[mask],{'pos':(3,)},comments)\ damask.Table(seeds[mask],{'pos':(3,)},comments)\
.add('material',materials[mask].astype(int))\ .add('material',material[mask].astype(int))\
.save(sys.stdout if name is None else os.path.splitext(name)[0]+'.seeds',legacy=True) .save(sys.stdout if name is None else os.path.splitext(name)[0]+'.seeds',legacy=True)

2
processing/pre/seeds_fromPokes.py
View File

@ -76,7 +76,7 @@ for name in filenames:
g[2] = k + offset[2] g[2] = k + offset[2]
g %= geom.grid g %= geom.grid
seeds[n,0:3] = (g+0.5)/geom.grid # normalize coordinates to box seeds[n,0:3] = (g+0.5)/geom.grid # normalize coordinates to box
seeds[n, 3] = geom.materials[g[0],g[1],g[2]] seeds[n, 3] = geom.material[g[0],g[1],g[2]]
if options.x: g[0] += 1 if options.x: g[0] += 1
if options.y: g[1] += 1 if options.y: g[1] += 1
n += 1 n += 1

250
python/damask/_geom.py
View File

@ -15,13 +15,13 @@ from . import grid_filters
class Geom: class Geom:
"""Geometry definition for grid solvers.""" """Geometry definition for grid solvers."""
def __init__(self,materials,size,origin=[0.0,0.0,0.0],comments=[]): def __init__(self,material,size,origin=[0.0,0.0,0.0],comments=[]):
""" """
New geometry definition from array of materials, size, and origin. New geometry definition from array of material, size, and origin.
Parameters Parameters
---------- ----------
materials : numpy.ndarray material : numpy.ndarray
Material index array (3D). Material index array (3D).
size : list or numpy.ndarray size : list or numpy.ndarray
Physical size of the geometry in meter. Physical size of the geometry in meter.
@ -31,16 +31,16 @@ class Geom:
Comment lines. Comment lines.
""" """
if len(materials.shape) != 3: if len(material.shape) != 3:
raise ValueError(f'Invalid materials shape {materials.shape}.') raise ValueError(f'Invalid material shape {material.shape}.')
elif materials.dtype not in np.sctypes['float'] + np.sctypes['int']: elif material.dtype not in np.sctypes['float'] + np.sctypes['int']:
raise TypeError(f'Invalid materials data type {materials.dtype}.') raise TypeError(f'Invalid material data type {material.dtype}.')
else: else:
self.materials = np.copy(materials) self.material = np.copy(material)
if self.materials.dtype in np.sctypes['float'] and \ if self.material.dtype in np.sctypes['float'] and \
np.all(self.materials == self.materials.astype(int).astype(float)): np.all(self.material == self.material.astype(int).astype(float)):
self.materials = self.materials.astype(int) self.material = self.material.astype(int)
if len(size) != 3 or any(np.array(size) <= 0): if len(size) != 3 or any(np.array(size) <= 0):
raise ValueError(f'Invalid size {size}.') raise ValueError(f'Invalid size {size}.')
@ -62,7 +62,7 @@ class Geom:
f'size x y z: {util.srepr(self.size, " x ")}', f'size x y z: {util.srepr(self.size, " x ")}',
f'origin x y z: {util.srepr(self.origin," ")}', f'origin x y z: {util.srepr(self.origin," ")}',
f'# materials: {self.N_materials}', f'# materials: {self.N_materials}',
f'max material: {np.nanmax(self.materials)}', f'max material: {np.nanmax(self.material)}',
]) ])
@ -103,21 +103,21 @@ class Geom:
message.append(util.delete(f'# materials: {other.N_materials}')) message.append(util.delete(f'# materials: {other.N_materials}'))
message.append(util.emph( f'# materials: { self.N_materials}')) message.append(util.emph( f'# materials: { self.N_materials}'))
if np.nanmax(other.materials) != np.nanmax(self.materials): if np.nanmax(other.material) != np.nanmax(self.material):
message.append(util.delete(f'max material: {np.nanmax(other.materials)}')) message.append(util.delete(f'max material: {np.nanmax(other.material)}'))
message.append(util.emph( f'max material: {np.nanmax( self.materials)}')) message.append(util.emph( f'max material: {np.nanmax( self.material)}'))
return util.return_message(message) return util.return_message(message)
@property @property
def grid(self): def grid(self):
return np.asarray(self.materials.shape) return np.asarray(self.material.shape)
@property @property
def N_materials(self): def N_materials(self):
return np.unique(self.materials).size return np.unique(self.material).size
@staticmethod @staticmethod
@ -160,7 +160,7 @@ class Geom:
else: else:
comments.append(line.strip()) comments.append(line.strip())
materials = np.empty(grid.prod()) # initialize as flat array material = np.empty(grid.prod()) # initialize as flat array
i = 0 i = 0
for line in content[header_length:]: for line in content[header_length:]:
items = line.split('#')[0].split() items = line.split('#')[0].split()
@ -172,16 +172,16 @@ class Geom:
abs(int(items[2])-int(items[0]))+1,dtype=float) abs(int(items[2])-int(items[0]))+1,dtype=float)
else: items = list(map(float,items)) else: items = list(map(float,items))
else: items = list(map(float,items)) else: items = list(map(float,items))
materials[i:i+len(items)] = items material[i:i+len(items)] = items
i += len(items) i += len(items)
if i != grid.prod(): if i != grid.prod():
raise TypeError(f'Invalid file: expected {grid.prod()} entries, found {i}') raise TypeError(f'Invalid file: expected {grid.prod()} entries, found {i}')
if not np.any(np.mod(materials,1) != 0.0): # no float present if not np.any(np.mod(material,1) != 0.0): # no float present
materials = materials.astype('int') material = material.astype('int')
return Geom(materials.reshape(grid,order='F'),size,origin,comments) return Geom(material.reshape(grid,order='F'),size,origin,comments)
@staticmethod @staticmethod
@ -200,9 +200,11 @@ class Geom:
comments = v.get_comments() comments = v.get_comments()
grid = np.array(v.vtk_data.GetDimensions())-1 grid = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
size = bbox[1] - bbox[0]
return Geom(v.get('material').reshape(grid,order='F'),size,bbox[0],comments=comments) return Geom(material = v.get('material').reshape(grid,order='F'),
size = bbox[1] - bbox[0],
origin = bbox[0],
comments=comments)
@staticmethod @staticmethod
@ -210,7 +212,7 @@ class Geom:
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights) return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
@staticmethod @staticmethod
def from_Laguerre_tessellation(grid,size,seeds,weights,materials=None,periodic=True): def from_Laguerre_tessellation(grid,size,seeds,weights,material=None,periodic=True):
""" """
Generate geometry from Laguerre tessellation. Generate geometry from Laguerre tessellation.
@ -224,7 +226,7 @@ class Geom:
Position of the seed points in meter. All points need to lay within the box. Position of the seed points in meter. All points need to lay within the box.
weights : numpy.ndarray of shape (seeds.shape[0]) weights : numpy.ndarray of shape (seeds.shape[0])
Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation. Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
materials : numpy.ndarray of shape (seeds.shape[0]), optional material : numpy.ndarray of shape (seeds.shape[0]), optional
Material ID of the seeds. Defaults to None, in which case materials are Material ID of the seeds. Defaults to None, in which case materials are
consecutively numbered. consecutively numbered.
periodic : Boolean, optional periodic : Boolean, optional
@ -246,27 +248,27 @@ class Geom:
result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords]) result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
pool.close() pool.close()
pool.join() pool.join()
materials_ = np.array(result.get()) material_ = np.array(result.get())
if periodic: if periodic:
materials_ = materials_.reshape(grid*3) material_ = material_.reshape(grid*3)
materials_ = materials_[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] material_ = material_[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
else: else:
materials_ = materials_.reshape(grid) material_ = material_.reshape(grid)
geom = Geom(materials = materials_+1, geom = Geom(material = material_+1,
size = size, size = size,
comments = util.execution_stamp('Geom','from_Laguerre_tessellation'), comments = util.execution_stamp('Geom','from_Laguerre_tessellation'),
) )
if materials is not None: if material is not None:
geom = geom.substitute(np.arange(seeds.shape[0])+1,materials) geom = geom.substitute(np.arange(seeds.shape[0])+1,material)
geom.comments = geom.comments[:-1] geom.comments = geom.comments[:-1]
return geom return geom
@staticmethod @staticmethod
def from_Voronoi_tessellation(grid,size,seeds,materials=None,periodic=True): def from_Voronoi_tessellation(grid,size,seeds,material=None,periodic=True):
""" """
Generate geometry from Voronoi tessellation. Generate geometry from Voronoi tessellation.
@ -278,7 +280,7 @@ class Geom:
Physical size of the geometry in meter. Physical size of the geometry in meter.
seeds : numpy.ndarray of shape (:,3) seeds : numpy.ndarray of shape (:,3)
Position of the seed points in meter. All points need to lay within the box. Position of the seed points in meter. All points need to lay within the box.
materials : numpy.ndarray of shape (seeds.shape[0]), optional material : numpy.ndarray of shape (seeds.shape[0]), optional
Material ID of the seeds. Defaults to None, in which case materials are Material ID of the seeds. Defaults to None, in which case materials are
consecutively numbered. consecutively numbered.
periodic : Boolean, optional periodic : Boolean, optional
@ -287,14 +289,14 @@ class Geom:
""" """
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3) coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,materials_ = KDTree.query(coords) devNull,material_ = KDTree.query(coords)
geom = Geom(materials = materials_.reshape(grid)+1, geom = Geom(material = material_.reshape(grid)+1,
size = size, size = size,
comments = util.execution_stamp('Geom','from_Voronoi_tessellation'), comments = util.execution_stamp('Geom','from_Voronoi_tessellation'),
) )
if materials is not None: if material is not None:
geom = geom.substitute(np.arange(seeds.shape[0])+1,materials) geom = geom.substitute(np.arange(seeds.shape[0])+1,material)
geom.comments = geom.comments[:-1] geom.comments = geom.comments[:-1]
return geom return geom
@ -327,10 +329,10 @@ class Geom:
plain = not compress plain = not compress
if plain: if plain:
format_string = '%g' if self.materials.dtype in np.sctypes['float'] else \ format_string = '%g' if self.material.dtype in np.sctypes['float'] else \
'%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.materials))))) '%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.material)))))
np.savetxt(fname, np.savetxt(fname,
self.materials.reshape([grid[0],np.prod(grid[1:])],order='F').T, self.material.reshape([grid[0],np.prod(grid[1:])],order='F').T,
header='\n'.join(header), fmt=format_string, comments='') header='\n'.join(header), fmt=format_string, comments='')
else: else:
try: try:
@ -341,7 +343,7 @@ class Geom:
compressType = None compressType = None
former = start = -1 former = start = -1
reps = 0 reps = 0
for current in self.materials.flatten('F'): for current in self.material.flatten('F'):
if abs(current - former) == 1 and (start - current) == reps*(former - current): if abs(current - former) == 1 and (start - current) == reps*(former - current):
compressType = 'to' compressType = 'to'
reps += 1 reps += 1
@ -386,7 +388,7 @@ class Geom:
""" """
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
v.add(self.materials.flatten(order='F'),'material') v.add(self.material.flatten(order='F'),'material')
v.add_comments(self.comments) v.add_comments(self.comments)
v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress) v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress)
@ -418,7 +420,7 @@ class Geom:
0 gives octahedron (|x|^(2^0) + |y|^(2^0) + |z|^(2^0) < 1) 0 gives octahedron (|x|^(2^0) + |y|^(2^0) + |z|^(2^0) < 1)
1 gives a sphere (|x|^(2^1) + |y|^(2^1) + |z|^(2^1) < 1) 1 gives a sphere (|x|^(2^1) + |y|^(2^1) + |z|^(2^1) < 1)
fill : int, optional fill : int, optional
Fill value for primitive. Defaults to materials.max() + 1. Fill value for primitive. Defaults to material.max() + 1.
R : damask.Rotation, optional R : damask.Rotation, optional
Rotation of primitive. Defaults to no rotation. Rotation of primitive. Defaults to no rotation.
inverse : Boolean, optional inverse : Boolean, optional
@ -444,12 +446,12 @@ class Geom:
if periodic: # translate back to center if periodic: # translate back to center
mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2)) mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))
fill_ = np.full_like(self.materials,np.nanmax(self.materials)+1 if fill is None else fill) fill_ = np.full_like(self.material,np.nanmax(self.material)+1 if fill is None else fill)
return Geom(materials = np.where(np.logical_not(mask) if inverse else mask, self.materials,fill_), return Geom(material = np.where(np.logical_not(mask) if inverse else mask, self.material,fill_),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','add_primitive')], comments = self.comments+[util.execution_stamp('Geom','add_primitive')],
) )
@ -471,19 +473,19 @@ class Geom:
raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.') raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
limits = [None,None] if reflect else [-2,0] limits = [None,None] if reflect else [-2,0]
ms = self.materials.copy() mat = self.material.copy()
if 'x' in directions: if 'x' in directions:
ms = np.concatenate([ms,ms[limits[0]:limits[1]:-1,:,:]],0) mat = np.concatenate([mat,mat[limits[0]:limits[1]:-1,:,:]],0)
if 'y' in directions: if 'y' in directions:
ms = np.concatenate([ms,ms[:,limits[0]:limits[1]:-1,:]],1) mat = np.concatenate([mat,mat[:,limits[0]:limits[1]:-1,:]],1)
if 'z' in directions: if 'z' in directions:
ms = np.concatenate([ms,ms[:,:,limits[0]:limits[1]:-1]],2) mat = np.concatenate([mat,mat[:,:,limits[0]:limits[1]:-1]],2)
return Geom(materials = ms, return Geom(material = mat,
size = self.size/self.grid*np.asarray(ms.shape), size = self.size/self.grid*np.asarray(mat.shape),
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','mirror')], comments = self.comments+[util.execution_stamp('Geom','mirror')],
) )
@ -502,12 +504,12 @@ class Geom:
if not set(directions).issubset(valid): if not set(directions).issubset(valid):
raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.') raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
ms = np.flip(self.materials, (valid.index(d) for d in directions if d in valid)) mat = np.flip(self.material, (valid.index(d) for d in directions if d in valid))
return Geom(materials = ms, return Geom(material = mat,
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','flip')], comments = self.comments+[util.execution_stamp('Geom','flip')],
) )
@ -523,17 +525,17 @@ class Geom:
Assume geometry to be periodic. Defaults to True. Assume geometry to be periodic. Defaults to True.
""" """
return Geom(materials = ndimage.interpolation.zoom( return Geom(material = ndimage.interpolation.zoom(
self.materials, self.material,
grid/self.grid, grid/self.grid,
output=self.materials.dtype, output=self.material.dtype,
order=0, order=0,
mode=('wrap' if periodic else 'nearest'), mode=('wrap' if periodic else 'nearest'),
prefilter=False prefilter=False
), ),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','scale')], comments = self.comments+[util.execution_stamp('Geom','scale')],
) )
@ -559,29 +561,29 @@ class Geom:
else: else:
return me return me
return Geom(materials = ndimage.filters.generic_filter( return Geom(material = ndimage.filters.generic_filter(
self.materials, self.material,
mostFrequent, mostFrequent,
size=(stencil if selection is None else stencil//2*2+1,)*3, size=(stencil if selection is None else stencil//2*2+1,)*3,
mode=('wrap' if periodic else 'nearest'), mode=('wrap' if periodic else 'nearest'),
extra_keywords=dict(selection=selection), extra_keywords=dict(selection=selection),
).astype(self.materials.dtype), ).astype(self.material.dtype),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','clean')], comments = self.comments+[util.execution_stamp('Geom','clean')],
) )
def renumber(self): def renumber(self):
"""Renumber sorted material indices to 1,...,N.""" """Renumber sorted material indices to 1,...,N."""
renumbered = np.empty(self.grid,dtype=self.materials.dtype) renumbered = np.empty(self.grid,dtype=self.material.dtype)
for i, oldID in enumerate(np.unique(self.materials)): for i, oldID in enumerate(np.unique(self.material)):
renumbered = np.where(self.materials == oldID, i+1, renumbered) renumbered = np.where(self.material == oldID, i+1, renumbered)
return Geom(materials = renumbered, return Geom(material = renumbered,
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','renumber')], comments = self.comments+[util.execution_stamp('Geom','renumber')],
) )
@ -594,32 +596,32 @@ class Geom:
R : damask.Rotation R : damask.Rotation
Rotation to apply to the geometry. Rotation to apply to the geometry.
fill : int or float, optional fill : int or float, optional
Material index to fill the corners. Defaults to materials.max() + 1. Material index to fill the corners. Defaults to material.max() + 1.
""" """
if fill is None: fill = np.nanmax(self.materials) + 1 if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if np.isnan(fill) or int(fill) != fill or self.materials.dtype==np.float else int dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int
Eulers = R.as_Eulers(degrees=True) Eulers = R.as_Eulers(degrees=True)
materials_in = self.materials.copy() material_in = self.material.copy()
# These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'') # These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
# see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf # see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
for angle,axes in zip(Eulers[::-1], [(0,1),(1,2),(0,1)]): for angle,axes in zip(Eulers[::-1], [(0,1),(1,2),(0,1)]):
materials_out = ndimage.rotate(materials_in,angle,axes,order=0, material_out = ndimage.rotate(material_in,angle,axes,order=0,
prefilter=False,output=dtype,cval=fill) prefilter=False,output=dtype,cval=fill)
if np.prod(materials_in.shape) == np.prod(materials_out.shape): if np.prod(material_in.shape) == np.prod(material_out.shape):
# avoid scipy interpolation errors for rotations close to multiples of 90° # avoid scipy interpolation errors for rotations close to multiples of 90°
materials_in = np.rot90(materials_in,k=np.rint(angle/90.).astype(int),axes=axes) material_in = np.rot90(material_in,k=np.rint(angle/90.).astype(int),axes=axes)
else: else:
materials_in = materials_out material_in = material_out
origin = self.origin-(np.asarray(materials_in.shape)-self.grid)*.5 * self.size/self.grid origin = self.origin-(np.asarray(material_in.shape)-self.grid)*.5 * self.size/self.grid
return Geom(materials = materials_in, return Geom(material = material_in,
size = self.size/self.grid*np.asarray(materials_in.shape), size = self.size/self.grid*np.asarray(material_in.shape),
origin = origin, origin = origin,
comments = self.comments+[util.execution_stamp('Geom','rotate')], comments = self.comments+[util.execution_stamp('Geom','rotate')],
) )
@ -634,12 +636,12 @@ class Geom:
offset : numpy.ndarray of shape (3) offset : numpy.ndarray of shape (3)
Offset (measured in grid points) from old to new geometry [0,0,0]. Offset (measured in grid points) from old to new geometry [0,0,0].
fill : int or float, optional fill : int or float, optional
Material index to fill the background. Defaults to materials.max() + 1. Material index to fill the background. Defaults to material.max() + 1.
""" """
if offset is None: offset = 0 if offset is None: offset = 0
if fill is None: fill = np.nanmax(self.materials) + 1 if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if int(fill) != fill or self.materials.dtype in np.sctypes['float'] else int dtype = float if int(fill) != fill or self.material.dtype in np.sctypes['float'] else int
canvas = np.full(self.grid if grid is None else grid,fill,dtype) canvas = np.full(self.grid if grid is None else grid,fill,dtype)
@ -648,35 +650,35 @@ class Geom:
ll = np.clip(-offset, 0,np.minimum( grid,self.grid-offset)) ll = np.clip(-offset, 0,np.minimum( grid,self.grid-offset))
ur = np.clip(-offset+self.grid,0,np.minimum( grid,self.grid-offset)) ur = np.clip(-offset+self.grid,0,np.minimum( grid,self.grid-offset))
canvas[ll[0]:ur[0],ll[1]:ur[1],ll[2]:ur[2]] = self.materials[LL[0]:UR[0],LL[1]:UR[1],LL[2]:UR[2]] canvas[ll[0]:ur[0],ll[1]:ur[1],ll[2]:ur[2]] = self.material[LL[0]:UR[0],LL[1]:UR[1],LL[2]:UR[2]]
return Geom(materials = canvas, return Geom(material = canvas,
size = self.size/self.grid*np.asarray(canvas.shape), size = self.size/self.grid*np.asarray(canvas.shape),
origin = self.origin+offset*self.size/self.grid, origin = self.origin+offset*self.size/self.grid,
comments = self.comments+[util.execution_stamp('Geom','canvas')], comments = self.comments+[util.execution_stamp('Geom','canvas')],
) )
def substitute(self,from_materials,to_materials): def substitute(self,from_material,to_material):
""" """
Substitute material indices. Substitute material indices.
Parameters Parameters
---------- ----------
from_materials : iterable of ints from_material : iterable of ints
Material indices to be substituted. Material indices to be substituted.
to_materials : iterable of ints to_material : iterable of ints
New material indices. New material indices.
""" """
substituted = self.materials.copy() substituted = self.material.copy()
for from_ms,to_ms in zip(from_materials,to_materials): for from_ms,to_ms in zip(from_material,to_material):
substituted[self.materials==from_ms] = to_ms substituted[self.material==from_ms] = to_ms
return Geom(materials = substituted, return Geom(material = substituted,
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','substitute')], comments = self.comments+[util.execution_stamp('Geom','substitute')],
) )
@ -695,7 +697,7 @@ class Geom:
Defaults to 1. Defaults to 1.
offset : int, optional offset : int, optional
Offset (positive or negative) to tag material indices, Offset (positive or negative) to tag material indices,
defaults to materials.max() + 1. defaults to material.max() + 1.
trigger : list of ints, optional trigger : list of ints, optional
List of material indices that trigger a change. List of material indices that trigger a change.
Defaults to [], meaning that any different neighbor triggers a change. Defaults to [], meaning that any different neighbor triggers a change.
@ -714,15 +716,15 @@ class Geom:
trigger = list(trigger) trigger = list(trigger)
return np.any(np.in1d(stencil,np.array(trigger))) return np.any(np.in1d(stencil,np.array(trigger)))
offset_ = np.nanmax(self.materials) if offset is None else offset offset_ = np.nanmax(self.material) if offset is None else offset
mask = ndimage.filters.generic_filter(self.materials, mask = ndimage.filters.generic_filter(self.material,
tainted_neighborhood, tainted_neighborhood,
size=1+2*vicinity, size=1+2*vicinity,
mode='wrap' if periodic else 'nearest', mode='wrap' if periodic else 'nearest',
extra_keywords={'trigger':trigger}) extra_keywords={'trigger':trigger})
return Geom(materials = np.where(mask, self.materials + offset_,self.materials), return Geom(material = np.where(mask, self.material + offset_,self.material),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','vicinity_offset')], comments = self.comments+[util.execution_stamp('Geom','vicinity_offset')],
) )

5
python/damask/util.py
View File

@ -117,6 +117,7 @@ def execute(cmd,
initialPath = os.getcwd() initialPath = os.getcwd()
myEnv = os.environ if env is None else env myEnv = os.environ if env is None else env
os.chdir(wd) os.chdir(wd)
print(f"executing '{cmd}' in '{wd}'")
process = subprocess.Popen(shlex.split(cmd), process = subprocess.Popen(shlex.split(cmd),
stdout = subprocess.PIPE, stdout = subprocess.PIPE,
stderr = subprocess.PIPE, stderr = subprocess.PIPE,
@ -128,7 +129,7 @@ def execute(cmd,
stdout = stdout.decode('utf-8').replace('\x08','') stdout = stdout.decode('utf-8').replace('\x08','')
stderr = stderr.decode('utf-8').replace('\x08','') stderr = stderr.decode('utf-8').replace('\x08','')
if process.returncode != 0: if process.returncode != 0:
raise RuntimeError(f'{cmd} failed with returncode {process.returncode}') raise RuntimeError(f"'{cmd}' failed with returncode {process.returncode}")
return stdout, stderr return stdout, stderr
@ -172,7 +173,7 @@ def scale_to_coprime(v):
m = (np.array(v) * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v)) ** 0.5).astype(np.int) m = (np.array(v) * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v)) ** 0.5).astype(np.int)
m = m//reduce(np.gcd,m) m = m//reduce(np.gcd,m)
with np.errstate(divide='ignore'): with np.errstate(invalid='ignore'):
if not np.allclose(np.ma.masked_invalid(v/m),v[np.argmax(abs(v))]/m[np.argmax(abs(v))]): if not np.allclose(np.ma.masked_invalid(v/m),v[np.argmax(abs(v))]/m[np.argmax(abs(v))]):
raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?') raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?')

52
python/tests/test_Geom.py
View File

@ -11,8 +11,8 @@ from damask import util
def geom_equal(a,b): def geom_equal(a,b):
return np.all(a.materials == b.materials) and \ return np.all(a.material == b.material) and \
np.all(a.grid == b.grid) and \ np.all(a.grid == b.grid) and \
np.allclose(a.size, b.size) and \ np.allclose(a.size, b.size) and \
str(a.diff(b)) == str(b.diff(a)) str(a.diff(b)) == str(b.diff(a))
@ -38,7 +38,7 @@ class TestGeom:
def test_diff_not_equal(self,default): def test_diff_not_equal(self,default):
new = Geom(default.materials[1:,1:,1:]+1,default.size*.9,np.ones(3)-default.origin,comments=['modified']) new = Geom(default.material[1:,1:,1:]+1,default.size*.9,np.ones(3)-default.origin,comments=['modified'])
assert str(default.diff(new)) != '' assert str(default.diff(new)) != ''
@ -93,28 +93,28 @@ class TestGeom:
def test_invalid_size(self,default): def test_invalid_size(self,default):
with pytest.raises(ValueError): with pytest.raises(ValueError):
Geom(default.materials[1:,1:,1:], Geom(default.material[1:,1:,1:],
size=np.ones(2)) size=np.ones(2))
def test_invalid_origin(self,default): def test_invalid_origin(self,default):
with pytest.raises(ValueError): with pytest.raises(ValueError):
Geom(default.materials[1:,1:,1:], Geom(default.material[1:,1:,1:],
size=np.ones(3), size=np.ones(3),
origin=np.ones(4)) origin=np.ones(4))
def test_invalid_materials_shape(self,default): def test_invalid_materials_shape(self,default):
materials = np.ones((3,3)) material = np.ones((3,3))
with pytest.raises(ValueError): with pytest.raises(ValueError):
Geom(materials, Geom(material,
size=np.ones(3)) size=np.ones(3))
def test_invalid_materials_type(self,default): def test_invalid_materials_type(self,default):
materials = np.random.randint(1,300,(3,4,5))==1 material = np.random.randint(1,300,(3,4,5))==1
with pytest.raises(TypeError): with pytest.raises(TypeError):
Geom(materials) Geom(material)
@pytest.mark.parametrize('directions,reflect',[ @pytest.mark.parametrize('directions,reflect',[
@ -205,10 +205,10 @@ class TestGeom:
def test_renumber(self,default): def test_renumber(self,default):
materials = default.materials.copy() material = default.material.copy()
for m in np.unique(materials): for m in np.unique(material):
materials[materials==m] = materials.max() + np.random.randint(1,30) material[material==m] = material.max() + np.random.randint(1,30)
modified = Geom(materials, modified = Geom(material,
default.size, default.size,
default.origin) default.origin)
assert not geom_equal(modified,default) assert not geom_equal(modified,default)
@ -218,13 +218,13 @@ class TestGeom:
def test_substitute(self,default): def test_substitute(self,default):
offset = np.random.randint(1,500) offset = np.random.randint(1,500)
modified = Geom(default.materials + offset, modified = Geom(default.material + offset,
default.size, default.size,
default.origin) default.origin)
assert not geom_equal(modified,default) assert not geom_equal(modified,default)
assert geom_equal(default, assert geom_equal(default,
modified.substitute(np.arange(default.materials.max())+1+offset, modified.substitute(np.arange(default.material.max())+1+offset,
np.arange(default.materials.max())+1)) np.arange(default.material.max())+1))
@pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]), @pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]),
@ -251,7 +251,7 @@ class TestGeom:
grid = default.grid grid = default.grid
grid_add = np.random.randint(0,30,(3)) grid_add = np.random.randint(0,30,(3))
modified = default.canvas(grid + grid_add) modified = default.canvas(grid + grid_add)
assert np.all(modified.materials[:grid[0],:grid[1],:grid[2]] == default.materials) assert np.all(modified.material[:grid[0],:grid[1],:grid[2]] == default.material)
@pytest.mark.parametrize('center1,center2',[(np.random.random(3)*.5,np.random.random()*8), @pytest.mark.parametrize('center1,center2',[(np.random.random(3)*.5,np.random.random()*8),
@ -271,7 +271,7 @@ class TestGeom:
s = np.random.random(3)+.5 s = np.random.random(3)+.5
G_1 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center1,exponent) G_1 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center1,exponent)
G_2 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center2,exponent) G_2 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center2,exponent)
assert np.count_nonzero(G_1.materials!=2) == np.count_nonzero(G_2.materials!=2) assert np.count_nonzero(G_1.material!=2) == np.count_nonzero(G_2.material!=2)
@pytest.mark.parametrize('center',[np.random.randint(4,10,(3)), @pytest.mark.parametrize('center',[np.random.randint(4,10,(3)),
@ -308,14 +308,14 @@ class TestGeom:
geom = Geom(m,np.random.rand(3)).vicinity_offset(vicinity,offset,trigger=trigger) geom = Geom(m,np.random.rand(3)).vicinity_offset(vicinity,offset,trigger=trigger)
assert np.all(m2==geom.materials) assert np.all(m2==geom.material)
@pytest.mark.parametrize('periodic',[True,False]) @pytest.mark.parametrize('periodic',[True,False])
def test_vicinity_offset_invariant(self,default,periodic): def test_vicinity_offset_invariant(self,default,periodic):
offset = default.vicinity_offset(trigger=[default.materials.max()+1, offset = default.vicinity_offset(trigger=[default.material.max()+1,
default.materials.min()-1]) default.material.min()-1])
assert np.all(offset.materials==default.materials) assert np.all(offset.material==default.material)
@pytest.mark.parametrize('periodic',[True,False]) @pytest.mark.parametrize('periodic',[True,False])
@ -338,7 +338,7 @@ class TestGeom:
ms = np.random.randint(1, N_seeds+1) ms = np.random.randint(1, N_seeds+1)
weights[ms-1] = np.random.random() weights[ms-1] = np.random.random()
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,periodic=np.random.random()>0.5) Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,periodic=np.random.random()>0.5)
assert np.all(Laguerre.materials == ms) assert np.all(Laguerre.material == ms)
@pytest.mark.parametrize('approach',['Laguerre','Voronoi']) @pytest.mark.parametrize('approach',['Laguerre','Voronoi'])
@ -346,10 +346,10 @@ class TestGeom:
grid = np.random.randint(5,10,3)*2 grid = np.random.randint(5,10,3)*2
size = grid.astype(np.float) size = grid.astype(np.float)
seeds = np.vstack((size*np.array([0.5,0.25,0.5]),size*np.array([0.5,0.75,0.5]))) seeds = np.vstack((size*np.array([0.5,0.25,0.5]),size*np.array([0.5,0.75,0.5])))
materials = np.ones(grid) material = np.ones(grid)
materials[:,grid[1]//2:,:] = 2 material[:,grid[1]//2:,:] = 2
if approach == 'Laguerre': if approach == 'Laguerre':
geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),periodic=np.random.random()>0.5) geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),periodic=np.random.random()>0.5)
elif approach == 'Voronoi': elif approach == 'Voronoi':
geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5) geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5)
assert np.all(geom.materials == materials) assert np.all(geom.material == material)

6
src/CPFEM.f90
View File

@ -106,7 +106,7 @@ subroutine CPFEM_init
num_commercialFEM, & num_commercialFEM, &
debug_CPFEM debug_CPFEM
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(6) print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
allocate(CPFEM_cs( 6,discretization_nIP,discretization_nElem), source= 0.0_pReal) allocate(CPFEM_cs( 6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
allocate(CPFEM_dcsdE( 6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal) allocate(CPFEM_dcsdE( 6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
@ -132,7 +132,7 @@ subroutine CPFEM_init
print'(a32,1x,6(i8,1x))', 'CPFEM_cs: ', shape(CPFEM_cs) print'(a32,1x,6(i8,1x))', 'CPFEM_cs: ', shape(CPFEM_cs)
print'(a32,1x,6(i8,1x))', 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE) print'(a32,1x,6(i8,1x))', 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
print'(a32,1x,6(i8,1x),/)', 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood) print'(a32,1x,6(i8,1x),/)', 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
flush(6) flush(IO_STDOUT)
endif endif
end subroutine CPFEM_init end subroutine CPFEM_init
@ -250,7 +250,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
'<< CPFEM >> stress/MPa at elFE ip ', elFE, ip, CPFEM_cs(1:6,ip,elCP)*1.0e-6_pReal '<< CPFEM >> stress/MPa at elFE ip ', elFE, ip, CPFEM_cs(1:6,ip,elCP)*1.0e-6_pReal
print'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))', & print'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))', &
'<< CPFEM >> Jacobian/GPa at elFE ip ', elFE, ip, transpose(CPFEM_dcsdE(1:6,1:6,ip,elCP))*1.0e-9_pReal '<< CPFEM >> Jacobian/GPa at elFE ip ', elFE, ip, transpose(CPFEM_dcsdE(1:6,1:6,ip,elCP))*1.0e-9_pReal
flush(6) flush(IO_STDOUT)
endif endif
endif endif

2
src/CPFEM2.f90
View File

@ -76,7 +76,7 @@ end subroutine CPFEM_initAll
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_init subroutine CPFEM_init
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(6) print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
if (interface_restartInc > 0) call crystallite_restartRead if (interface_restartInc > 0) call crystallite_restartRead

16
src/DAMASK_interface.f90
View File

@ -14,7 +14,7 @@
#define PETSC_MINOR_MAX 13 #define PETSC_MINOR_MAX 13
module DAMASK_interface module DAMASK_interface
use, intrinsic :: iso_fortran_env use, intrinsic :: ISO_fortran_env
use PETScSys use PETScSys
@ -82,7 +82,7 @@ subroutine DAMASK_interface_init
print'(/,a)', ' <<<+- DAMASK_interface init -+>>>' print'(/,a)', ' <<<+- DAMASK_interface init -+>>>'
open(6, encoding='UTF-8') ! for special characters in output open(OUTPUT_unit, encoding='UTF-8') ! for special characters in output
! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK%203 ! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK%203
#ifdef DEBUG #ifdef DEBUG
@ -101,8 +101,8 @@ subroutine DAMASK_interface_init
#endif #endif
print*, achar(27)//'[0m' print*, achar(27)//'[0m'
print'(a)', ' Roters et al., Computational Materials Science 158:420–478, 2019' print*, 'Roters et al., Computational Materials Science 158:420–478, 2019'
print'(a)', ' https://doi.org/10.1016/j.commatsci.2018.04.030' print*, 'https://doi.org/10.1016/j.commatsci.2018.04.030'
print'(/,a)', ' Version: '//DAMASKVERSION print'(/,a)', ' Version: '//DAMASKVERSION
@ -373,7 +373,7 @@ function makeRelativePath(a,b)
a_cleaned = rectifyPath(trim(a)//'/') a_cleaned = rectifyPath(trim(a)//'/')
b_cleaned = rectifyPath(b) b_cleaned = rectifyPath(b)
do i = 1, min(1024,len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned))) do i = 1, min(len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned)))
if (a_cleaned(i:i) /= b_cleaned(i:i)) exit if (a_cleaned(i:i) /= b_cleaned(i:i)) exit
if (a_cleaned(i:i) == '/') posLastCommonSlash = i if (a_cleaned(i:i) == '/') posLastCommonSlash = i
enddo enddo
@ -395,7 +395,7 @@ subroutine catchSIGTERM(signal) bind(C)
integer(C_INT), value :: signal integer(C_INT), value :: signal
interface_SIGTERM = .true. interface_SIGTERM = .true.
print'(a,i2.2,a)', ' received signal ',signal, ', set SIGTERM=TRUE' print'(a,i0,a)', ' received signal ',signal, ', set SIGTERM=TRUE'
end subroutine catchSIGTERM end subroutine catchSIGTERM
@ -420,7 +420,7 @@ subroutine catchSIGUSR1(signal) bind(C)
integer(C_INT), value :: signal integer(C_INT), value :: signal
interface_SIGUSR1 = .true. interface_SIGUSR1 = .true.
print'(a,i2.2,a)', ' received signal ',signal, ', set SIGUSR1=TRUE' print'(a,i0,a)', ' received signal ',signal, ', set SIGUSR1=TRUE'
end subroutine catchSIGUSR1 end subroutine catchSIGUSR1
@ -445,7 +445,7 @@ subroutine catchSIGUSR2(signal) bind(C)
integer(C_INT), value :: signal integer(C_INT), value :: signal
interface_SIGUSR2 = .true. interface_SIGUSR2 = .true.
print'(a,i2.2,a)', ' received signal ',signal, ', set SIGUSR2=TRUE' print'(a,i0,a)', ' received signal ',signal, ', set SIGUSR2=TRUE'
end subroutine catchSIGUSR2 end subroutine catchSIGUSR2

2
src/DAMASK_marc.f90
View File

@ -30,7 +30,7 @@
module DAMASK_interface module DAMASK_interface
use prec use prec
#if __INTEL_COMPILER >= 1800 #if __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: & use, intrinsic :: ISO_fortran_env, only: &
compiler_version, & compiler_version, &
compiler_options compiler_options
#endif #endif

61
src/IO.f90
View File

@ -6,6 +6,10 @@
!> @brief input/output functions !> @brief input/output functions
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module IO module IO
use, intrinsic :: ISO_fortran_env, only: &
IO_STDOUT => OUTPUT_UNIT, &
IO_STDERR => ERROR_UNIT
use prec use prec
implicit none implicit none
@ -16,7 +20,7 @@ module IO
character, parameter, public :: & character, parameter, public :: &
IO_EOL = new_line('DAMASK'), & !< end of line character IO_EOL = new_line('DAMASK'), & !< end of line character
IO_COMMENT = '#' IO_COMMENT = '#'
character(len=*), parameter, private :: & character(len=*), parameter :: &
IO_DIVIDER = '───────────────────'//& IO_DIVIDER = '───────────────────'//&
'───────────────────'//& '───────────────────'//&
'───────────────────'//& '───────────────────'//&
@ -37,7 +41,8 @@ module IO
IO_stringAsFloat, & IO_stringAsFloat, &
IO_stringAsBool, & IO_stringAsBool, &
IO_error, & IO_error, &
IO_warning IO_warning, &
IO_STDOUT
contains contains
@ -47,7 +52,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine IO_init subroutine IO_init
print'(/,a)', ' <<<+- IO init -+>>>'; flush(6) print'(/,a)', ' <<<+- IO init -+>>>'; flush(IO_STDOUT)
call selfTest call selfTest
@ -538,29 +543,29 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
end select end select
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(0,'(/,a)') ' ┌'//IO_DIVIDER//'┐' write(IO_STDERR,'(/,a)') ' ┌'//IO_DIVIDER//'┐'
write(0,'(a,24x,a,40x,a)') ' │','error', '│' write(IO_STDERR,'(a,24x,a,40x,a)') ' │','error', '│'
write(0,'(a,24x,i3,42x,a)') ' │',error_ID, '│' write(IO_STDERR,'(a,24x,i3,42x,a)') ' │',error_ID, '│'
write(0,'(a)') ' ├'//IO_DIVIDER//'┤' write(IO_STDERR,'(a)') ' ├'//IO_DIVIDER//'┤'
write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(msg)),',',& write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(msg)),',',&
max(1,72-len_trim(msg)-4),'x,a)' max(1,72-len_trim(msg)-4),'x,a)'
write(0,formatString) '│ ',trim(msg), '│' write(IO_STDERR,formatString) '│ ',trim(msg), '│'
if (present(ext_msg)) then if (present(ext_msg)) then
write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',& write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
max(1,72-len_trim(ext_msg)-4),'x,a)' max(1,72-len_trim(ext_msg)-4),'x,a)'
write(0,formatString) '│ ',trim(ext_msg), '│' write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│'
endif endif
if (present(el)) & if (present(el)) &
write(0,'(a19,1x,i9,44x,a3)') ' │ at element ',el, '│' write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at element ',el, '│'
if (present(ip)) & if (present(ip)) &
write(0,'(a19,1x,i9,44x,a3)') ' │ at IP ',ip, '│' write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at IP ',ip, '│'
if (present(g)) & if (present(g)) &
write(0,'(a19,1x,i9,44x,a3)') ' │ at constituent',g, '│' write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at constituent',g, '│'
if (present(instance)) & if (present(instance)) &
write(0,'(a19,1x,i9,44x,a3)') ' │ at instance ',instance, '│' write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at instance ',instance, '│'
write(0,'(a,69x,a)') ' │', '│' write(IO_STDERR,'(a,69x,a)') ' │', '│'
write(0,'(a)') ' └'//IO_DIVIDER//'┘' write(IO_STDERR,'(a)') ' └'//IO_DIVIDER//'┘'
flush(0) flush(IO_STDERR)
call quit(9000+error_ID) call quit(9000+error_ID)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -623,27 +628,27 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
end select end select
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(/,a)') ' ┌'//IO_DIVIDER//'┐' write(IO_STDERR,'(/,a)') ' ┌'//IO_DIVIDER//'┐'
write(6,'(a,24x,a,38x,a)') ' │','warning', '│' write(IO_STDERR,'(a,24x,a,38x,a)') ' │','warning', '│'
write(6,'(a,24x,i3,42x,a)') ' │',warning_ID, '│' write(IO_STDERR,'(a,24x,i3,42x,a)') ' │',warning_ID, '│'
write(6,'(a)') ' ├'//IO_DIVIDER//'┤' write(IO_STDERR,'(a)') ' ├'//IO_DIVIDER//'┤'
write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(msg)),',',& write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(msg)),',',&
max(1,72-len_trim(msg)-4),'x,a)' max(1,72-len_trim(msg)-4),'x,a)'
write(6,formatString) '│ ',trim(msg), '│' write(IO_STDERR,formatString) '│ ',trim(msg), '│'
if (present(ext_msg)) then if (present(ext_msg)) then
write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',& write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
max(1,72-len_trim(ext_msg)-4),'x,a)' max(1,72-len_trim(ext_msg)-4),'x,a)'
write(6,formatString) '│ ',trim(ext_msg), '│' write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│'
endif endif
if (present(el)) & if (present(el)) &
write(6,'(a19,1x,i9,44x,a3)') ' │ at element ',el, '│' write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at element ',el, '│'
if (present(ip)) & if (present(ip)) &
write(6,'(a19,1x,i9,44x,a3)') ' │ at IP ',ip, '│' write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at IP ',ip, '│'
if (present(g)) & if (present(g)) &
write(6,'(a19,1x,i9,44x,a3)') ' │ at constituent',g, '│' write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at constituent',g, '│'
write(6,'(a,69x,a)') ' │', '│' write(IO_STDERR,'(a,69x,a)') ' │', '│'
write(6,'(a)') ' └'//IO_DIVIDER//'┘' write(IO_STDERR,'(a)') ' └'//IO_DIVIDER//'┘'
flush(6) flush(IO_STDERR)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
end subroutine IO_warning end subroutine IO_warning

2
src/base64.f90
View File

@ -27,7 +27,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine base64_init subroutine base64_init
print'(/,a)', ' <<<+- base64 init -+>>>'; flush(6) print'(/,a)', ' <<<+- base64 init -+>>>'; flush(IO_STDOUT)
call selfTest call selfTest

8
src/config.f90
View File

@ -35,7 +35,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine config_init subroutine config_init
print'(/,a)', ' <<<+- config init -+>>>'; flush(6) print'(/,a)', ' <<<+- config init -+>>>'; flush(IO_STDOUT)
call parse_material call parse_material
call parse_numerics call parse_numerics
@ -59,7 +59,7 @@ subroutine parse_material
inquire(file=fname,exist=fileExists) inquire(file=fname,exist=fileExists)
if(.not. fileExists) call IO_error(100,ext_msg=fname) if(.not. fileExists) call IO_error(100,ext_msg=fname)
endif endif
print*, 'reading '//fname; flush(6) print*, 'reading '//fname; flush(IO_STDOUT)
config_material => YAML_parse_file(fname) config_material => YAML_parse_file(fname)
end subroutine parse_material end subroutine parse_material
@ -75,7 +75,7 @@ subroutine parse_numerics
config_numerics => emptyDict config_numerics => emptyDict
inquire(file='numerics.yaml', exist=fexist) inquire(file='numerics.yaml', exist=fexist)
if (fexist) then if (fexist) then
print*, 'reading numerics.yaml'; flush(6) print*, 'reading numerics.yaml'; flush(IO_STDOUT)
config_numerics => YAML_parse_file('numerics.yaml') config_numerics => YAML_parse_file('numerics.yaml')
endif endif
@ -92,7 +92,7 @@ subroutine parse_debug
config_debug => emptyDict config_debug => emptyDict
inquire(file='debug.yaml', exist=fexist) inquire(file='debug.yaml', exist=fexist)
fileExists: if (fexist) then fileExists: if (fexist) then
print*, 'reading debug.yaml'; flush(6) print*, 'reading debug.yaml'; flush(IO_STDOUT)
config_debug => YAML_parse_file('debug.yaml') config_debug => YAML_parse_file('debug.yaml')
endif fileExists endif fileExists

2
src/constitutive.f90
View File

@ -446,7 +446,7 @@ subroutine constitutive_init
call damage_init call damage_init
call thermal_init call thermal_init
print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(6) print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT)
constitutive_source_maxSizeDotState = 0 constitutive_source_maxSizeDotState = 0
PhaseLoop2:do p = 1,phases%length PhaseLoop2:do p = 1,phases%length

2
src/constitutive_plastic_disloTungsten.f90
View File

@ -100,7 +100,7 @@ module function plastic_disloTungsten_init() result(myPlasticity)
myPlasticity = plastic_active('disloTungsten') myPlasticity = plastic_active('disloTungsten')
Ninstance = count(myPlasticity) Ninstance = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016'

2
src/constitutive_plastic_dislotwin.f90
View File

@ -147,7 +147,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
myPlasticity = plastic_active('dislotwin') myPlasticity = plastic_active('dislotwin')
Ninstance = count(myPlasticity) Ninstance = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'

2
src/constitutive_plastic_isotropic.f90
View File

@ -71,7 +71,7 @@ module function plastic_isotropic_init() result(myPlasticity)
myPlasticity = plastic_active('isotropic') myPlasticity = plastic_active('isotropic')
Ninstance = count(myPlasticity) Ninstance = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'

2
src/constitutive_plastic_kinehardening.f90
View File

@ -83,7 +83,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
myPlasticity = plastic_active('kinehardening') myPlasticity = plastic_active('kinehardening')
Ninstance = count(myPlasticity) Ninstance = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
allocate(param(Ninstance)) allocate(param(Ninstance))

2
src/constitutive_plastic_none.f90
View File

@ -35,7 +35,7 @@ module function plastic_none_init() result(myPlasticity)
enddo enddo
Ninstance = count(myPlasticity) Ninstance = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
do p = 1, phases%length do p = 1, phases%length

28
src/constitutive_plastic_nonlocal.f90
View File

@ -189,7 +189,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
myPlasticity = plastic_active('nonlocal') myPlasticity = plastic_active('nonlocal')
Ninstance = count(myPlasticity) Ninstance = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) then if(Ninstance == 0) then
call geometry_plastic_nonlocal_disable call geometry_plastic_nonlocal_disable
return return
@ -199,7 +199,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
print*, 'Kords, Dissertation RWTH Aachen, 2014' print*, 'Kords, Dissertation RWTH Aachen, 2014'
print*, 'http://publications.rwth-aachen.de/record/229993'//IO_EOL print*, 'http://publications.rwth-aachen.de/record/229993'
allocate(param(Ninstance)) allocate(param(Ninstance))
allocate(state(Ninstance)) allocate(state(Ninstance))
@ -741,10 +741,10 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
if (debugConstitutive%extensive & if (debugConstitutive%extensive &
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)& .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
.or. .not. debugConstitutive%selective)) then .or. .not. debugConstitutive%selective)) then
write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_microstructure at el ip ',el,ip print'(/,a,i8,1x,i2,1x,i1,/)', '<< CONST >> nonlocal_microstructure at el ip ',el,ip
write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', stt%rho_forest(:,of) print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,of)
write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,of)*1e-6 print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,of)*1e-6
write(6,'(a,/,12x,12(f10.5,1x),/)') '<< CONST >> tauBack / MPa', dst%tau_back(:,of)*1e-6 print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,of)*1e-6
endif endif
#endif #endif
@ -958,8 +958,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
if (debugConstitutive%extensive & if (debugConstitutive%extensive &
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)& .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
.or. .not. debugConstitutive%selective)) then .or. .not. debugConstitutive%selective)) then
write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation remobilization', deltaRhoRemobilization(:,1:8) print'(a,/,8(12x,12(e12.5,1x),/))', '<< CONST >> dislocation remobilization', deltaRhoRemobilization(:,1:8)
write(6,'(a,/,10(12x,12(e12.5,1x),/),/)') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress print'(a,/,10(12x,12(e12.5,1x),/),/)', '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress
endif endif
#endif #endif
@ -1047,8 +1047,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
if (debugConstitutive%basic & if (debugConstitutive%basic &
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip) & .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip) &
.or. .not. debugConstitutive%selective)) then .or. .not. debugConstitutive%selective)) then
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip print'(a,/,10(12x,12(e12.5,1x),/))', '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot print'(a,/,4(12x,12(e12.5,1x),/))', '<< CONST >> gdot / 1/s',gdot
endif endif
#endif #endif
@ -1156,8 +1156,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
.or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then
#ifdef DEBUG #ifdef DEBUG
if (debugConstitutive%extensive) then if (debugConstitutive%extensive) then
write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
write(6,'(a)') '<< CONST >> enforcing cutback !!!' print'(a)', '<< CONST >> enforcing cutback !!!'
endif endif
#endif #endif
plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
@ -1268,8 +1268,8 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG #ifdef DEBUG
if (debugConstitutive%extensive) then if (debugConstitutive%extensive) then
write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(gdot) > 0.0_pReal & maxval(abs(v0), abs(gdot) > 0.0_pReal &
.and. prm%CFLfactor * abs(v0) * timestep & .and. prm%CFLfactor * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &

2
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -92,7 +92,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
myPlasticity = plastic_active('phenopowerlaw') myPlasticity = plastic_active('phenopowerlaw')
Ninstance = count(myPlasticity) Ninstance = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
allocate(param(Ninstance)) allocate(param(Ninstance))

4
src/crystallite.f90
View File

@ -294,7 +294,7 @@ subroutine crystallite_init
print'(a42,1x,i10)', ' # of elements: ', eMax print'(a42,1x,i10)', ' # of elements: ', eMax
print'(a42,1x,i10)', ' # of integration points/element: ', iMax print'(a42,1x,i10)', ' # of integration points/element: ', iMax
print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif
@ -1561,7 +1561,7 @@ subroutine crystallite_restartWrite
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: fileName, datasetName character(len=pStringLen) :: fileName, datasetName
print*, ' writing field and constitutive data required for restart to file';flush(6) print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT)
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'a') fileHandle = HDF5_openFile(fileName,'a')

2
src/damage_local.f90
View File

@ -49,7 +49,7 @@ subroutine damage_local_init
homog, & homog, &
homogDamage homogDamage
print'(/,a)', ' <<<+- damage_local init -+>>>'; flush(6) print'(/,a)', ' <<<+- damage_local init -+>>>'; flush(IO_STDOUT)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
! read numerics parameter and do sanity check ! read numerics parameter and do sanity check

2
src/element.f90
View File

@ -922,7 +922,7 @@ subroutine tElement_init(self,elemType)
self%nIPneighbors = size(self%IPneighbor,1) self%nIPneighbors = size(self%IPneighbor,1)
print'(/,a)', ' <<<+- element_init -+>>>'; flush(6) print'(/,a)', ' <<<+- element_init -+>>>'; flush(IO_STDOUT)
print*, 'element type: ',self%elemType print*, 'element type: ',self%elemType
print*, ' geom type: ',self%geomType print*, ' geom type: ',self%geomType

67
src/grid/DAMASK_grid.f90
View File

@ -99,10 +99,10 @@ program DAMASK_grid
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll call CPFEM_initAll
write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>'; flush(6) print'(/,a)', ' <<<+- DAMASK_spectral init -+>>>'; flush(IO_STDOUT)
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize field solver information ! initialize field solver information
@ -263,56 +263,56 @@ program DAMASK_grid
reportAndCheck: if (worldrank == 0) then reportAndCheck: if (worldrank == 0) then
write (loadcase_string, '(i0)' ) currentLoadCase write (loadcase_string, '(i0)' ) currentLoadCase
write(6,'(/,1x,a,i0)') 'load case: ', currentLoadCase print'(/,a,i0)', ' load case: ', currentLoadCase
if (.not. newLoadCase%followFormerTrajectory) & if (.not. newLoadCase%followFormerTrajectory) &
write(6,'(2x,a)') 'drop guessing along trajectory' print*, ' drop guessing along trajectory'
if (newLoadCase%deformation%myType == 'l') then if (newLoadCase%deformation%myType == 'l') then
do j = 1, 3 do j = 1, 3
if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. &
any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined
enddo enddo
write(6,'(2x,a)') 'velocity gradient:' print*, ' velocity gradient:'
else if (newLoadCase%deformation%myType == 'f') then else if (newLoadCase%deformation%myType == 'f') then
write(6,'(2x,a)') 'deformation gradient at end of load case:' print*, ' deformation gradient at end of load case:'
else else
write(6,'(2x,a)') 'deformation gradient rate:' print*, ' deformation gradient rate:'
endif endif
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if(newLoadCase%deformation%maskLogical(i,j)) then if(newLoadCase%deformation%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j) write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j)
else else
write(6,'(2x,12a)',advance='no') ' * ' write(IO_STDOUT,'(2x,12a)',advance='no') ' * '
endif endif
enddo; write(6,'(/)',advance='no') enddo; write(IO_STDOUT,'(/)',advance='no')
enddo enddo
if (any(newLoadCase%stress%maskLogical .eqv. & if (any(newLoadCase%stress%maskLogical .eqv. &
newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only
if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) & if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) &
.and. (math_I3<1))) errorID = 838 ! no rotation is allowed by stress BC .and. (math_I3<1))) errorID = 838 ! no rotation is allowed by stress BC
write(6,'(2x,a)') 'stress / GPa:' print*, ' stress / GPa:'
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if(newLoadCase%stress%maskLogical(i,j)) then if(newLoadCase%stress%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal
else else
write(6,'(2x,12a)',advance='no') ' * ' write(IO_STDOUT,'(2x,12a)',advance='no') ' * '
endif endif
enddo; write(6,'(/)',advance='no') enddo; write(IO_STDOUT,'(/)',advance='no')
enddo enddo
if (any(abs(matmul(newLoadCase%rot%asMatrix(), & if (any(abs(matmul(newLoadCase%rot%asMatrix(), &
transpose(newLoadCase%rot%asMatrix()))-math_I3) > & transpose(newLoadCase%rot%asMatrix()))-math_I3) > &
reshape(spread(tol_math_check,1,9),[ 3,3]))) errorID = 846 ! given rotation matrix contains strain reshape(spread(tol_math_check,1,9),[ 3,3]))) errorID = 846 ! given rotation matrix contains strain
if (any(dNeq(newLoadCase%rot%asMatrix(), math_I3))) & if (any(dNeq(newLoadCase%rot%asMatrix(), math_I3))) &
write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
transpose(newLoadCase%rot%asMatrix()) transpose(newLoadCase%rot%asMatrix())
if (newLoadCase%time < 0.0_pReal) errorID = 834 ! negative time increment if (newLoadCase%time < 0.0_pReal) errorID = 834 ! negative time increment
write(6,'(2x,a,f0.3)') 'time: ', newLoadCase%time print'(a,f0.3)', ' time: ', newLoadCase%time
if (newLoadCase%incs < 1) errorID = 835 ! non-positive incs count if (newLoadCase%incs < 1) errorID = 835 ! non-positive incs count
write(6,'(2x,a,i0)') 'increments: ', newLoadCase%incs print'(a,i0)', ' increments: ', newLoadCase%incs
if (newLoadCase%outputfrequency < 1) errorID = 836 ! non-positive result frequency if (newLoadCase%outputfrequency < 1) errorID = 836 ! non-positive result frequency
write(6,'(2x,a,i0)') 'output frequency: ', newLoadCase%outputfrequency print'(a,i0)', ' output frequency: ', newLoadCase%outputfrequency
if (newLoadCase%restartfrequency < 1) errorID = 839 ! non-positive restart frequency if (newLoadCase%restartfrequency < 1) errorID = 839 ! non-positive restart frequency
if (newLoadCase%restartfrequency < huge(0)) & if (newLoadCase%restartfrequency < huge(0)) &
write(6,'(2x,a,i0)') 'restart frequency: ', newLoadCase%restartfrequency print'(a,i0)', ' restart frequency: ', newLoadCase%restartfrequency
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
endif reportAndCheck endif reportAndCheck
loadCases = [loadCases,newLoadCase] ! load case is ok, append it loadCases = [loadCases,newLoadCase] ! load case is ok, append it
@ -341,9 +341,8 @@ program DAMASK_grid
writeHeader: if (interface_restartInc < 1) then writeHeader: if (interface_restartInc < 1) then
open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
if (debug_grid%contains('basic')) & if (debug_grid%contains('basic')) print'(/,a)', ' header of statistics file written out'
write(6,'(/,a)') ' header of statistics file written out' flush(IO_STDOUT)
flush(6)
else writeHeader else writeHeader
open(newunit=statUnit,file=trim(getSolverJobName())//& open(newunit=statUnit,file=trim(getSolverJobName())//&
'.sta',form='FORMATTED', position='APPEND', status='OLD') '.sta',form='FORMATTED', position='APPEND', status='OLD')
@ -351,7 +350,7 @@ program DAMASK_grid
endif endif
writeUndeformed: if (interface_restartInc < 1) then writeUndeformed: if (interface_restartInc < 1) then
write(6,'(1/,a)') ' ... writing initial configuration to file ........................' print'(/,a)', ' ... writing initial configuration to file ........................'
call CPFEM_results(0,0.0_pReal) call CPFEM_results(0,0.0_pReal)
endif writeUndeformed endif writeUndeformed
@ -397,8 +396,8 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
write(6,'(/,a)') ' ###########################################################################' print'(/,a)', ' ###########################################################################'
write(6,'(1x,a,es12.5,6(a,i0))') & print'(1x,a,es12.5,6(a,i0))', &
'Time', time, & 'Time', time, &
's: Increment ', inc,'/',loadCases(currentLoadCase)%incs,& 's: Increment ', inc,'/',loadCases(currentLoadCase)%incs,&
'-', stepFraction,'/',subStepFactor**cutBackLevel,& '-', stepFraction,'/',subStepFactor**cutBackLevel,&
@ -406,7 +405,7 @@ program DAMASK_grid
write(incInfo,'(4(a,i0))') & write(incInfo,'(4(a,i0))') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-', stepFraction,'/',subStepFactor**cutBackLevel '-', stepFraction,'/',subStepFactor**cutBackLevel
flush(6) flush(IO_STDOUT)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forward fields ! forward fields
@ -475,7 +474,7 @@ program DAMASK_grid
cutBackLevel = cutBackLevel + 1 cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time time = time - timeinc ! rewind time
timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep
write(6,'(/,a)') ' cutting back ' print'(/,a)', ' cutting back '
else ! no more options to continue else ! no more options to continue
call IO_warning(850) call IO_warning(850)
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
@ -487,14 +486,14 @@ program DAMASK_grid
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then if (all(solres(:)%converged)) then
write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' converged' print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged'
else else
write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' NOT converged' print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6) endif; flush(IO_STDOUT)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................' print'(1/,a)', ' ... writing results to file ......................................'
flush(6) flush(IO_STDOUT)
call CPFEM_results(totalIncsCounter,time) call CPFEM_results(totalIncsCounter,time)
endif endif
if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then
@ -510,7 +509,7 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report summary of whole calculation ! report summary of whole calculation
write(6,'(/,a)') ' ###########################################################################' print'(/,a)', ' ###########################################################################'
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
call quit(0) ! no complains ;) call quit(0) ! no complains ;)

54
src/grid/discretization_grid.f90
View File

@ -56,7 +56,7 @@ subroutine discretization_grid_init(restart)
myGrid !< domain grid of this process myGrid !< domain grid of this process
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
microstructureAt materialAt
integer :: & integer :: &
j, & j, &
@ -65,12 +65,12 @@ subroutine discretization_grid_init(restart)
integer(C_INTPTR_T) :: & integer(C_INTPTR_T) :: &
devNull, z, z_offset devNull, z, z_offset
write(6,'(/,a)') ' <<<+- discretization_grid init -+>>>'; flush(6) print'(/,a)', ' <<<+- discretization_grid init -+>>>'; flush(IO_STDOUT)
if(index(interface_geomFile,'.vtr') /= 0) then if(index(interface_geomFile,'.vtr') /= 0) then
call readVTR(grid,geomSize,origin,microstructureAt) call readVTR(grid,geomSize,origin,materialAt)
else else
call readGeom(grid,geomSize,origin,microstructureAt) call readGeom(grid,geomSize,origin,materialAt)
endif endif
print'(/,a,3(i12 ))', ' grid a b c: ', grid print'(/,a,3(i12 ))', ' grid a b c: ', grid
@ -102,15 +102,14 @@ subroutine discretization_grid_init(restart)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! general discretization ! general discretization
microstructureAt = microstructureAt(product(grid(1:2))*grid3Offset+1: & materialAt = materialAt(product(grid(1:2))*grid3Offset+1:product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI
product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI
call discretization_init(microstructureAt, & call discretization_init(materialAt, &
IPcoordinates0(myGrid,mySize,grid3Offset), & IPcoordinates0(myGrid,mySize,grid3Offset), &
Nodes0(myGrid,mySize,grid3Offset),& Nodes0(myGrid,mySize,grid3Offset),&
merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write bottom layer merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write top layer...
(grid(1)+1) * (grid(2)+1) * grid3,& ! do not write bottom layer (is top of rank-1) (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process
worldrank<1)) worldrank+1==worldsize))
FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements
FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP
@ -147,7 +146,7 @@ end subroutine discretization_grid_init
!> @details important variables have an implicit "save" attribute. Therefore, this function is !> @details important variables have an implicit "save" attribute. Therefore, this function is
! supposed to be called only once! ! supposed to be called only once!
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine readGeom(grid,geomSize,origin,microstructure) subroutine readGeom(grid,geomSize,origin,material)
integer, dimension(3), intent(out) :: & integer, dimension(3), intent(out) :: &
grid ! grid (across all processes!) grid ! grid (across all processes!)
@ -155,7 +154,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure)
geomSize, & ! size (across all processes!) geomSize, & ! size (across all processes!)
origin ! origin (across all processes!) origin ! origin (across all processes!)
integer, dimension(:), intent(out), allocatable :: & integer, dimension(:), intent(out), allocatable :: &
microstructure material
character(len=:), allocatable :: rawData character(len=:), allocatable :: rawData
character(len=65536) :: line character(len=65536) :: line
@ -167,7 +166,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure)
startPos, endPos, & startPos, endPos, &
myStat, & myStat, &
l, & !< line counter l, & !< line counter
c, & !< counter for # microstructures in line c, & !< counter for # materials in line
o, & !< order of "to" packing o, & !< order of "to" packing
e, & !< "element", i.e. spectral collocation point e, & !< "element", i.e. spectral collocation point
i, j i, j
@ -266,7 +265,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure)
if(any(geomSize < 0.0_pReal)) & if(any(geomSize < 0.0_pReal)) &
call IO_error(error_ID = 842, ext_msg='size (readGeom)') call IO_error(error_ID = 842, ext_msg='size (readGeom)')
allocate(microstructure(product(grid)), source = -1) ! too large in case of MPI (shrink later, not very elegant) allocate(material(product(grid)), source = -1) ! too large in case of MPI (shrink later, not very elegant)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! read and interpret content ! read and interpret content
@ -281,18 +280,18 @@ subroutine readGeom(grid,geomSize,origin,microstructure)
noCompression: if (chunkPos(1) /= 3) then noCompression: if (chunkPos(1) /= 3) then
c = chunkPos(1) c = chunkPos(1)
microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)] material(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)]
else noCompression else noCompression
compression: if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'of') then compression: if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'of') then
c = IO_intValue(line,chunkPos,1) c = IO_intValue(line,chunkPos,1)
microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,3),i = 1,IO_intValue(line,chunkPos,1))] material(e:e+c-1) = [(IO_intValue(line,chunkPos,3),i = 1,IO_intValue(line,chunkPos,1))]
else if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to') then compression else if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to') then compression
c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1 c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1
o = merge(+1, -1, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1)) o = merge(+1, -1, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1))
microstructure(e:e+c-1) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)] material(e:e+c-1) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)]
else compression else compression
c = chunkPos(1) c = chunkPos(1)
microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)] material(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)]
endif compression endif compression
endif noCompression endif noCompression
@ -308,7 +307,7 @@ end subroutine readGeom
!> @brief Parse vtk rectilinear grid (.vtr) !> @brief Parse vtk rectilinear grid (.vtr)
!> @details https://vtk.org/Wiki/VTK_XML_Formats !> @details https://vtk.org/Wiki/VTK_XML_Formats
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine readVTR(grid,geomSize,origin,microstructure) subroutine readVTR(grid,geomSize,origin,material)
integer, dimension(3), intent(out) :: & integer, dimension(3), intent(out) :: &
grid ! grid (across all processes!) grid ! grid (across all processes!)
@ -316,7 +315,7 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
geomSize, & ! size (across all processes!) geomSize, & ! size (across all processes!)
origin ! origin (across all processes!) origin ! origin (across all processes!)
integer, dimension(:), intent(out), allocatable :: & integer, dimension(:), intent(out), allocatable :: &
microstructure material
character(len=:), allocatable :: fileContent, dataType, headerType character(len=:), allocatable :: fileContent, dataType, headerType
logical :: inFile,inGrid,gotCoordinates,gotCellData,compressed logical :: inFile,inGrid,gotCoordinates,gotCellData,compressed
@ -364,11 +363,9 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
else else
if(index(fileContent(startPos:endPos),'<CellData>',kind=pI64) /= 0_pI64) then if(index(fileContent(startPos:endPos),'<CellData>',kind=pI64) /= 0_pI64) then
gotCellData = .true. gotCellData = .true.
startPos = endPos + 2_pI64
do while (index(fileContent(startPos:endPos),'</CellData>',kind=pI64) == 0_pI64) do while (index(fileContent(startPos:endPos),'</CellData>',kind=pI64) == 0_pI64)
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
if(index(fileContent(startPos:endPos),'<DataArray',kind=pI64) /= 0_pI64 .and. & if(index(fileContent(startPos:endPos),'<DataArray',kind=pI64) /= 0_pI64 .and. &
getXMLValue(fileContent(startPos:endPos),'Name') == 'materialpoint' ) then getXMLValue(fileContent(startPos:endPos),'Name') == 'material' ) then
if(getXMLValue(fileContent(startPos:endPos),'format') /= 'binary') & if(getXMLValue(fileContent(startPos:endPos),'format') /= 'binary') &
call IO_error(error_ID = 844, ext_msg='format (materialpoint)') call IO_error(error_ID = 844, ext_msg='format (materialpoint)')
@ -377,10 +374,11 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
startPos = endPos + 2_pI64 startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64 endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64 ! start (no leading whitespace) s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64 ! start (no leading whitespace)
microstructure = as_Int(fileContent(s:endPos),headerType,compressed,dataType) material = as_Int(fileContent(s:endPos),headerType,compressed,dataType)
exit exit
endif endif
startPos = endPos + 2_pI64 startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
enddo enddo
elseif(index(fileContent(startPos:endPos),'<Coordinates>',kind=pI64) /= 0_pI64) then elseif(index(fileContent(startPos:endPos),'<Coordinates>',kind=pI64) /= 0_pI64) then
gotCoordinates = .true. gotCoordinates = .true.
@ -415,10 +413,10 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
end do end do
if(.not. allocated(microstructure)) call IO_error(error_ID = 844, ext_msg='materialpoint not found') if(.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found')
if(size(microstructure) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(materialpoint)') if(size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)')
if(any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size') if(any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if(any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid') if(any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid')
contains contains

42
src/grid/grid_damage_spectral.f90
View File

@ -22,42 +22,38 @@ module grid_damage_spectral
implicit none implicit none
private private
type, private :: tNumerics type :: tNumerics
integer :: & integer :: &
itmax !< max number of iterations itmax !< maximum number of iterations
real(pReal) :: & real(pReal) :: &
residualStiffness, & !< non-zero residual damage residualStiffness, & !< non-zero residual damage
eps_damage_atol, & !< absolute tolerance for damage evolution eps_damage_atol, & !< absolute tolerance for damage evolution
eps_damage_rtol !< relative tolerance for damage evolution eps_damage_rtol !< relative tolerance for damage evolution
end type tNumerics end type tNumerics
type(tNumerics), private :: num type(tNumerics) :: num
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams), private :: params
type(tSolutionParams) :: params
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
SNES, private :: damage_snes SNES :: damage_snes
Vec, private :: solution_vec Vec :: solution_vec
PetscInt, private :: xstart, xend, ystart, yend, zstart, zend PetscInt :: xstart, xend, ystart, yend, zstart, zend
real(pReal), private, dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
phi_current, & !< field of current damage phi_current, & !< field of current damage
phi_lastInc, & !< field of previous damage phi_lastInc, & !< field of previous damage
phi_stagInc !< field of staggered damage phi_stagInc !< field of staggered damage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc. ! reference diffusion tensor, mobility etc.
integer, private :: totalIter = 0 !< total iteration in current increment integer :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3), private :: K_ref real(pReal), dimension(3,3) :: K_ref
real(pReal), private :: mu_ref real(pReal) :: mu_ref
public :: & public :: &
grid_damage_spectral_init, & grid_damage_spectral_init, &
grid_damage_spectral_solution, & grid_damage_spectral_solution, &
grid_damage_spectral_forward grid_damage_spectral_forward
private :: &
formResidual
contains contains
@ -77,10 +73,10 @@ subroutine grid_damage_spectral_init
character(len=pStringLen) :: & character(len=pStringLen) :: &
snes_type snes_type
write(6,'(/,a)') ' <<<+- grid_spectral_damage init -+>>>' print'(/,a)', ' <<<+- grid_spectral_damage init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
@ -152,8 +148,6 @@ subroutine grid_damage_spectral_init
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr) call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! damage reference diffusion update
call updateReference call updateReference
end subroutine grid_damage_spectral_init end subroutine grid_damage_spectral_init
@ -210,11 +204,11 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr) call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr)
if (solution%converged) & if (solution%converged) &
write(6,'(/,a)') ' ... nonlocal damage converged .....................................' print'(/,a)', ' ... nonlocal damage converged .....................................'
write(6,'(/,a,f8.6,2x,f8.6,2x,e11.4,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& write(IO_STDOUT,'(/,a,f8.6,2x,f8.6,2x,e11.4,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',&
phi_min, phi_max, stagNorm phi_min, phi_max, stagNorm
write(6,'(/,a)') ' ===========================================================================' print'(/,a)', ' ==========================================================================='
flush(6) flush(IO_STDOUT)
end function grid_damage_spectral_solution end function grid_damage_spectral_solution

37
src/grid/grid_mech_FEM.f90
View File

@ -122,7 +122,7 @@ subroutine grid_mech_FEM_init
PetscScalar, pointer, dimension(:,:,:,:) :: & PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc u_current,u_lastInc
write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>'; flush(6) print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT)
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! debugging options ! debugging options
@ -130,13 +130,12 @@ subroutine grid_mech_FEM_init
debugRotation = debug_grid%contains('rotation') debugRotation = debug_grid%contains('rotation')
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameter and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict) num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal) num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) num%itmin = num_grid%get_asInt ('itmin',defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
@ -225,7 +224,7 @@ subroutine grid_mech_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
restartRead: if (interface_restartInc > 0) then restartRead: if (interface_restartInc > 0) then
write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file' print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file'
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName) fileHandle = HDF5_openFile(fileName)
@ -254,7 +253,7 @@ subroutine grid_mech_FEM_init
CHKERRQ(ierr) CHKERRQ(ierr)
restartRead2: if (interface_restartInc > 0) then restartRead2: if (interface_restartInc > 0) then
write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file' print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
call HDF5_read(groupHandle,C_volAvg, 'C_volAvg') call HDF5_read(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
@ -304,11 +303,11 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr);CHKERRQ(ierr) call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check convergence ! check convergence
call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr) call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr)
solution%converged = reason > 0 solution%converged = reason > 0
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
@ -353,7 +352,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!-------------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = & F_aimDot = &
@ -414,7 +413,7 @@ subroutine grid_mech_FEM_restartWrite
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
write(6,'(a)') ' writing solver data required for restart to file'; flush(6) print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'w') fileHandle = HDF5_openFile(fileName,'w')
@ -476,13 +475,13 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
write(6,'(1/,a)') ' ... reporting .............................................................' print'(1/,a)', ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')'
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & print'(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ===========================================================================' print'(/,a)', ' ==========================================================================='
flush(6) flush(IO_STDOUT)
end subroutine converged end subroutine converged
@ -516,13 +515,13 @@ subroutine formResidual(da_local,x_local, &
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax
if (debugRotation) & if (debugRotation) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim) ' deformation gradient aim =', transpose(F_aim)
flush(6) flush(IO_STDOUT)
endif newIteration endif newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -541,7 +540,7 @@ subroutine formResidual(da_local,x_local, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call Utilities_constitutiveResponse(P_current,& call utilities_constitutiveResponse(P_current,&
P_av,C_volAvg,devNull, & P_av,C_volAvg,devNull, &
F,params%timeinc,params%rotation_BC) F,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)

38
src/grid/grid_mech_spectral_basic.f90
View File

@ -42,8 +42,7 @@ module grid_mech_spectral_basic
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
logical, private:: & logical, private :: debugRotation
debugRotation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
@ -110,13 +109,13 @@ subroutine grid_mech_spectral_basic_init
character(len=pStringLen) :: & character(len=pStringLen) :: &
fileName fileName
write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(6) print'(/,a)', ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(IO_STDOUT)
write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' print*, 'Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012' print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! debugging options ! debugging options
@ -132,7 +131,6 @@ subroutine grid_mech_spectral_basic_init
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.01_pReal) num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.01_pReal)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) num%itmin = num_grid%get_asInt ('itmin',defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
@ -186,7 +184,7 @@ subroutine grid_mech_spectral_basic_init
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! places pointer on PETSc data call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
restartRead: if (interface_restartInc > 0) then restartRead: if (interface_restartInc > 0) then
write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file' print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file'
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName) fileHandle = HDF5_openFile(fileName)
@ -211,7 +209,7 @@ subroutine grid_mech_spectral_basic_init
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead2: if (interface_restartInc > 0) then restartRead2: if (interface_restartInc > 0) then
write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file' print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
call HDF5_read(groupHandle,C_volAvg, 'C_volAvg') call HDF5_read(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
@ -377,7 +375,7 @@ subroutine grid_mech_spectral_basic_restartWrite
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
write(6,'(a)') ' writing solver data required for restart to file'; flush(6) print'(a)', ' writing solver data required for restart to file'; flush(IO_STDOUT)
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'w') fileHandle = HDF5_openFile(fileName,'w')
@ -437,13 +435,13 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
write(6,'(1/,a)') ' ... reporting .............................................................' print'(1/,a)', ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')'
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & print'(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ===========================================================================' print'(/,a)', ' ==========================================================================='
flush(6) flush(IO_STDOUT)
end subroutine converged end subroutine converged
@ -475,13 +473,13 @@ subroutine formResidual(in, F, &
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if (debugRotation) & if (debugRotation) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim) ' deformation gradient aim =', transpose(F_aim)
flush(6) flush(IO_STDOUT)
endif newIteration endif newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

35
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -123,10 +123,10 @@ subroutine grid_mech_spectral_polarisation_init
character(len=pStringLen) :: & character(len=pStringLen) :: &
fileName fileName
write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(6) print'(/,a)', ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(IO_STDOUT)
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
! debugging options ! debugging options
@ -134,9 +134,8 @@ subroutine grid_mech_spectral_polarisation_init
debugRotation = debug_grid%contains('rotation') debugRotation = debug_grid%contains('rotation')
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters ! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict) num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
@ -207,7 +206,7 @@ subroutine grid_mech_spectral_polarisation_init
F_tau => FandF_tau(9:17,:,:,:) F_tau => FandF_tau(9:17,:,:,:)
restartRead: if (interface_restartInc > 0) then restartRead: if (interface_restartInc > 0) then
write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file' print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file'
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName) fileHandle = HDF5_openFile(fileName)
@ -236,7 +235,7 @@ subroutine grid_mech_spectral_polarisation_init
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead2: if (interface_restartInc > 0) then restartRead2: if (interface_restartInc > 0) then
write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file' print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
call HDF5_read(groupHandle,C_volAvg, 'C_volAvg') call HDF5_read(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
@ -434,7 +433,7 @@ subroutine grid_mech_spectral_polarisation_restartWrite
F => FandF_tau(0: 8,:,:,:) F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:) F_tau => FandF_tau(9:17,:,:,:)
write(6,'(a)') ' writing solver data required for restart to file'; flush(6) print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'w') fileHandle = HDF5_openFile(fileName,'w')
@ -498,15 +497,15 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
write(6,'(1/,a)') ' ... reporting .............................................................' print'(1/,a)', ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & print '(a,f12.2,a,es8.2,a,es9.2,a)', ' error curl = ', &
err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & print '(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ===========================================================================' print'(/,a)', ' ==========================================================================='
flush(6) flush(IO_STDOUT)
end subroutine converged end subroutine converged
@ -558,13 +557,13 @@ subroutine formResidual(in, FandF_tau, &
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if(debugRotation) & if(debugRotation) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim) ' deformation gradient aim =', transpose(F_aim)
flush(6) flush(IO_STDOUT)
endif newIteration endif newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

34
src/grid/grid_thermal_spectral.f90
View File

@ -23,20 +23,17 @@ module grid_thermal_spectral
implicit none implicit none
private private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams) :: params
type :: tNumerics type :: tNumerics
integer :: & integer :: &
itmax !< maximum number of iterations itmax !< maximum number of iterations
real(pReal) :: & real(pReal) :: &
eps_thermal_atol, & !< absolute tolerance for thermal equilibrium eps_thermal_atol, & !< absolute tolerance for thermal equilibrium
eps_thermal_rtol !< relative tolerance for thermal equilibrium eps_thermal_rtol !< relative tolerance for thermal equilibrium
end type tNumerics end type tNumerics
type(tNumerics) :: num type(tNumerics) :: num
type(tSolutionParams) :: params
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
SNES :: thermal_snes SNES :: thermal_snes
@ -74,13 +71,13 @@ subroutine grid_thermal_spectral_init
class(tNode), pointer :: & class(tNode), pointer :: &
num_grid num_grid
write(6,'(/,a)') ' <<<+- grid_thermal_spectral init -+>>>' print'(/,a)', ' <<<+- grid_thermal_spectral init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameter and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict) num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
num%eps_thermal_atol = num_grid%get_asFloat ('eps_thermal_atol',defaultVal=1.0e-2_pReal) num%eps_thermal_atol = num_grid%get_asFloat ('eps_thermal_atol',defaultVal=1.0e-2_pReal)
@ -94,8 +91,7 @@ subroutine grid_thermal_spectral_init
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,& call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
num_grid%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -110,7 +106,7 @@ subroutine grid_thermal_spectral_init
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid grid(1),grid(2),grid(3), & ! global grid
1, 1, worldsize, & 1, 1, worldsize, &
1, 0, & ! #dof (thermal phase field), ghost boundary width (domain overlap) 1, 0, & ! #dof (T field), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & ! local grid [grid(1)],[grid(2)],localK, & ! local grid
thermal_grid,ierr) ! handle, error thermal_grid,ierr) ! handle, error
CHKERRQ(ierr) CHKERRQ(ierr)
@ -159,8 +155,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
timeinc_old !< increment in time of last increment timeinc_old !< increment in time of last increment
integer :: i, j, k, cell integer :: i, j, k, cell
type(tSolutionState) :: solution type(tSolutionState) :: solution
class(tNode), pointer :: &
num_grid
PetscInt :: devNull PetscInt :: devNull
PetscReal :: T_min, T_max, stagNorm, solnNorm PetscReal :: T_min, T_max, stagNorm, solnNorm
@ -204,11 +198,11 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr) call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr)
if (solution%converged) & if (solution%converged) &
write(6,'(/,a)') ' ... thermal conduction converged ..................................' print'(/,a)', ' ... thermal conduction converged ..................................'
write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',& write(IO_STDOUT,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',&
T_min, T_max, stagNorm T_min, T_max, stagNorm
write(6,'(/,a)') ' ===========================================================================' print'(/,a)', ' ==========================================================================='
flush(6) flush(IO_STDOUT)
end function grid_thermal_spectral_solution end function grid_thermal_spectral_solution

44
src/grid/spectral_utilities.f90
View File

@ -208,10 +208,10 @@ subroutine spectral_utilities_init
debugPETSc = debug_grid%contains('petsc') debugPETSc = debug_grid%contains('petsc')
if(debugPETSc) write(6,'(3(/,a),/)') & if(debugPETSc) print'(3(/,a),/)', &
' Initializing PETSc with debug options: ', & ' Initializing PETSc with debug options: ', &
trim(PETScDebug), & trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.yaml '; flush(6) ' add more using the PETSc_Options keyword in numerics.yaml '; flush(IO_STDOUT)
num_grid => config_numerics%get('grid',defaultVal=emptyDict) num_grid => config_numerics%get('grid',defaultVal=emptyDict)
@ -280,7 +280,7 @@ subroutine spectral_utilities_init
if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match' if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match'
call fftw_set_timelimit(num_grid%get_asFloat('fftw_timelimit',defaultVal=-1.0_pReal)) call fftw_set_timelimit(num_grid%get_asFloat('fftw_timelimit',defaultVal=-1.0_pReal))
if (debugGeneral) write(6,'(/,a)') ' FFTW initialized'; flush(6) print*, 'FFTW initialized'; flush(IO_STDOUT)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! MPI allocation ! MPI allocation
@ -506,8 +506,8 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
logical :: err logical :: err
write(6,'(/,a)') ' ... doing gamma convolution ...............................................' print'(/,a)', ' ... doing gamma convolution ...............................................'
flush(6) flush(IO_STDOUT)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium) ! do the actual spectral method calculation (mechanical equilibrium)
@ -576,8 +576,8 @@ real(pReal) function utilities_divergenceRMS()
integer :: i, j, k, ierr integer :: i, j, k, ierr
complex(pReal), dimension(3) :: rescaledGeom complex(pReal), dimension(3) :: rescaledGeom
write(6,'(/,a)') ' ... calculating divergence ................................................' print'(/,a)', ' ... calculating divergence ................................................'
flush(6) flush(IO_STDOUT)
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
@ -620,8 +620,8 @@ real(pReal) function utilities_curlRMS()
complex(pReal), dimension(3,3) :: curl_fourier complex(pReal), dimension(3,3) :: curl_fourier
complex(pReal), dimension(3) :: rescaledGeom complex(pReal), dimension(3) :: rescaledGeom
write(6,'(/,a)') ' ... calculating curl ......................................................' print'(/,a)', ' ... calculating curl ......................................................'
flush(6) flush(IO_STDOUT)
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
@ -700,10 +700,10 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
temp99_real = math_3333to99(rot_BC%rotate(C)) temp99_real = math_3333to99(rot_BC%rotate(C))
if(debugGeneral) then if(debugGeneral) then
write(6,'(/,a)') ' ... updating masked compliance ............................................' print'(/,a)', ' ... updating masked compliance ............................................'
write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& write(IO_STDOUT,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',&
transpose(temp99_Real)*1.0e-9_pReal transpose(temp99_Real)*1.0e-9_pReal
flush(6) flush(IO_STDOUT)
endif endif
do i = 1,9; do j = 1,9 do i = 1,9; do j = 1,9
@ -723,9 +723,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
if (debugGeneral .or. errmatinv) then if (debugGeneral .or. errmatinv) then
write(formatString, '(i2)') size_reduced write(formatString, '(i2)') size_reduced
formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
write(6,trim(formatString),advance='no') ' C * S (load) ', & write(IO_STDOUT,trim(formatString),advance='no') ' C * S (load) ', &
transpose(matmul(c_reduced,s_reduced)) transpose(matmul(c_reduced,s_reduced))
write(6,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced) write(IO_STDOUT,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced)
if(errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance') if(errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
endif endif
temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9]) temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9])
@ -736,9 +736,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
utilities_maskedCompliance = math_99to3333(temp99_Real) utilities_maskedCompliance = math_99to3333(temp99_Real)
if(debugGeneral) then if(debugGeneral) then
write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') & write(IO_STDOUT,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
flush(6) flush(IO_STDOUT)
endif endif
end function utilities_maskedCompliance end function utilities_maskedCompliance
@ -822,8 +822,8 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
real(pReal) :: dPdF_norm_max, dPdF_norm_min real(pReal) :: dPdF_norm_max, dPdF_norm_min
real(pReal), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF real(pReal), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF
write(6,'(/,a)') ' ... evaluating constitutive response ......................................' print'(/,a)', ' ... evaluating constitutive response ......................................'
flush(6) flush(IO_STDOUT)
materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
@ -833,13 +833,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (debugRotation) & if (debugRotation) &
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
transpose(P_av)*1.e-6_pReal transpose(P_av)*1.e-6_pReal
if(present(rotation_BC)) & if(present(rotation_BC)) &
P_av = rotation_BC%rotate(P_av) P_av = rotation_BC%rotate(P_av)
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
transpose(P_av)*1.e-6_pReal transpose(P_av)*1.e-6_pReal
flush(6) flush(IO_STDOUT)
dPdF_max = 0.0_pReal dPdF_max = 0.0_pReal
dPdF_norm_max = 0.0_pReal dPdF_norm_max = 0.0_pReal
@ -1095,7 +1095,7 @@ subroutine utilities_saveReferenceStiffness
fileUnit,ierr fileUnit,ierr
if (worldrank == 0) then if (worldrank == 0) then
write(6,'(a)') ' writing reference stiffness data required for restart to file'; flush(6) print'(a)', ' writing reference stiffness data required for restart to file'; flush(IO_STDOUT)
open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',& open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',&
status='replace',access='stream',action='write',iostat=ierr) status='replace',access='stream',action='write',iostat=ierr)
if(ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref') if(ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref')

2
src/homogenization.f90
View File

@ -186,7 +186,7 @@ subroutine homogenization_init
materialpoint_F = materialpoint_F0 ! initialize to identity materialpoint_F = materialpoint_F0 ! initialize to identity
allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal)
print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(6) print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT)
num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10) num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10)
num%subStepMinHomog = num_homogGeneric%get_asFloat('subStepMin', defaultVal=1.0e-3_pReal) num%subStepMinHomog = num_homogGeneric%get_asFloat('subStepMin', defaultVal=1.0e-3_pReal)

28
src/homogenization_mech_RGC.f90
View File

@ -95,7 +95,7 @@ module subroutine mech_RGC_init(num_homogMech)
print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>' print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID) Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
print*, 'Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' print*, 'Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
@ -247,7 +247,7 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
print'(1x,3(e15.8,1x))',(F(i,j,iGrain), j = 1,3) print'(1x,3(e15.8,1x))',(F(i,j,iGrain), j = 1,3)
enddo enddo
print*,' ' print*,' '
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif
enddo enddo
@ -376,7 +376,7 @@ module procedure mech_RGC_updateState
'@ grain ',stresLoc(3),' in component ',stresLoc(1),stresLoc(2) '@ grain ',stresLoc(3),' in component ',stresLoc(1),stresLoc(2)
print'(a,e15.8,a,i3,a,i2)',' Max residual: ',residMax, & print'(a,e15.8,a,i3,a,i2)',' Max residual: ',residMax, &
' @ iface ',residLoc(1),' in direction ',residLoc(2) ' @ iface ',residLoc(1),' in direction ',residLoc(2)
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif
@ -388,7 +388,7 @@ module procedure mech_RGC_updateState
mech_RGC_updateState = .true. mech_RGC_updateState = .true.
#ifdef DEBUG #ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) & if (debugHomog%extensive .and. prm%of_debug == of) &
print*, '... done and happy'; flush(6) print*, '... done and happy'; flush(IO_STDOUT)
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -416,7 +416,7 @@ module procedure mech_RGC_updateState
print'(a,e15.8,/)', ' Volume discrepancy: ', dst%volumeDiscrepancy(of) print'(a,e15.8,/)', ' Volume discrepancy: ', dst%volumeDiscrepancy(of)
print'(a,e15.8)', ' Maximum relaxation rate: ', dst%relaxationRate_max(of) print'(a,e15.8)', ' Maximum relaxation rate: ', dst%relaxationRate_max(of)
print'(a,e15.8,/)', ' Average relaxation rate: ', dst%relaxationRate_avg(of) print'(a,e15.8,/)', ' Average relaxation rate: ', dst%relaxationRate_avg(of)
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif
@ -429,7 +429,7 @@ module procedure mech_RGC_updateState
#ifdef DEBUG #ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) & if (debugHomog%extensive .and. prm%of_debug == of) &
print'(a,/)', ' ... broken'; flush(6) print'(a,/)', ' ... broken'; flush(IO_STDOUT)
#endif #endif
return return
@ -437,7 +437,7 @@ module procedure mech_RGC_updateState
else ! proceed with computing the Jacobian and state update else ! proceed with computing the Jacobian and state update
#ifdef DEBUG #ifdef DEBUG
if (debugHomog%extensive .and. prm%of_debug == of) & if (debugHomog%extensive .and. prm%of_debug == of) &
print'(a,/)', ' ... not yet done'; flush(6) print'(a,/)', ' ... not yet done'; flush(IO_STDOUT)
#endif #endif
endif endif
@ -499,7 +499,7 @@ module procedure mech_RGC_updateState
print'(1x,100(e11.4,1x))',(smatrix(i,j), j = 1,3*nIntFaceTot) print'(1x,100(e11.4,1x))',(smatrix(i,j), j = 1,3*nIntFaceTot)
enddo enddo
print*,' ' print*,' '
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif
@ -559,7 +559,7 @@ module procedure mech_RGC_updateState
print'(1x,100(e11.4,1x))',(pmatrix(i,j), j = 1,3*nIntFaceTot) print'(1x,100(e11.4,1x))',(pmatrix(i,j), j = 1,3*nIntFaceTot)
enddo enddo
print*,' ' print*,' '
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif
@ -578,7 +578,7 @@ module procedure mech_RGC_updateState
print'(1x,100(e11.4,1x))',(rmatrix(i,j), j = 1,3*nIntFaceTot) print'(1x,100(e11.4,1x))',(rmatrix(i,j), j = 1,3*nIntFaceTot)
enddo enddo
print*,' ' print*,' '
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif
@ -593,7 +593,7 @@ module procedure mech_RGC_updateState
print'(1x,100(e11.4,1x))',(jmatrix(i,j), j = 1,3*nIntFaceTot) print'(1x,100(e11.4,1x))',(jmatrix(i,j), j = 1,3*nIntFaceTot)
enddo enddo
print*,' ' print*,' '
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif
@ -609,7 +609,7 @@ module procedure mech_RGC_updateState
print'(1x,100(e11.4,1x))',(jnverse(i,j), j = 1,3*nIntFaceTot) print'(1x,100(e11.4,1x))',(jnverse(i,j), j = 1,3*nIntFaceTot)
enddo enddo
print*,' ' print*,' '
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif
@ -625,7 +625,7 @@ module procedure mech_RGC_updateState
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback' print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback'
print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax)) print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
flush(6) flush(IO_STDOUT)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -636,7 +636,7 @@ module procedure mech_RGC_updateState
print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of) print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of)
enddo enddo
print*,' ' print*,' '
flush(6) flush(IO_STDOUT)
endif endif
#endif #endif

2
src/homogenization_mech_isostrain.f90
View File

@ -40,7 +40,7 @@ module subroutine mech_isostrain_init
print'(/,a)', ' <<<+- homogenization_mech_isostrain init -+>>>' print'(/,a)', ' <<<+- homogenization_mech_isostrain init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) Ninstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
allocate(param(Ninstance)) ! one container of parameters per instance allocate(param(Ninstance)) ! one container of parameters per instance

2
src/homogenization_mech_none.f90
View File

@ -21,7 +21,7 @@ module subroutine mech_none_init
print'(/,a)', ' <<<+- homogenization_mech_none init -+>>>' print'(/,a)', ' <<<+- homogenization_mech_none init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID) Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
do h = 1, size(homogenization_type) do h = 1, size(homogenization_type)
if (homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle if (homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle

2
src/kinematics_cleavage_opening.f90
View File

@ -49,7 +49,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
myKinematics = kinematics_active('cleavage_opening',kinematics_length) myKinematics = kinematics_active('cleavage_opening',kinematics_length)
Ninstance = count(myKinematics) Ninstance = count(myKinematics)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')

2
src/kinematics_slipplane_opening.f90
View File

@ -52,7 +52,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
myKinematics = kinematics_active('slipplane_opening',kinematics_length) myKinematics = kinematics_active('slipplane_opening',kinematics_length)
Ninstance = count(myKinematics) Ninstance = count(myKinematics)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')

2
src/kinematics_thermal_expansion.f90
View File

@ -42,7 +42,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
myKinematics = kinematics_active('thermal_expansion',kinematics_length) myKinematics = kinematics_active('thermal_expansion',kinematics_length)
Ninstance = count(myKinematics) Ninstance = count(myKinematics)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')

2
src/lattice.f90
View File

@ -457,7 +457,7 @@ subroutine lattice_init
phase, & phase, &
elasticity elasticity
print'(/,a)', ' <<<+- lattice init -+>>>'; flush(6) print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT)
phases => config_material%get('phase') phases => config_material%get('phase')
Nphases = phases%length Nphases = phases%length

41
src/marc/discretization_marc.f90
View File

@ -52,7 +52,7 @@ subroutine discretization_marc_init
type(tElement) :: elem type(tElement) :: elem
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
microstructureAt materialAt
integer:: & integer:: &
Nnodes, & !< total number of nodes in the mesh Nnodes, & !< total number of nodes in the mesh
Nelems, & !< total number of elements in the mesh Nelems, & !< total number of elements in the mesh
@ -70,7 +70,7 @@ subroutine discretization_marc_init
class(tNode), pointer :: & class(tNode), pointer :: &
num_commercialFEM num_commercialFEM
write(6,'(/,a)') ' <<<+- discretization_marc init -+>>>'; flush(6) print'(/,a)', ' <<<+- discretization_marc init -+>>>'; flush(6)
!--------------------------------------------------------------------------------- !---------------------------------------------------------------------------------
! read debug parameters ! read debug parameters
@ -83,7 +83,7 @@ subroutine discretization_marc_init
mesh_unitlength = num_commercialFEM%get_asFloat('unitlength',defaultVal=1.0_pReal) ! set physical extent of a length unit in mesh mesh_unitlength = num_commercialFEM%get_asFloat('unitlength',defaultVal=1.0_pReal) ! set physical extent of a length unit in mesh
if (mesh_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength') if (mesh_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength')
call inputRead(elem,node0_elem,connectivity_elem,microstructureAt) call inputRead(elem,node0_elem,connectivity_elem,materialAt)
nElems = size(connectivity_elem,2) nElems = size(connectivity_elem,2)
if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element')
@ -103,7 +103,7 @@ subroutine discretization_marc_init
call buildIPcoordinates(IP_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,& call buildIPcoordinates(IP_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,&
elem%nIPs*nElems]),node0_cell) elem%nIPs*nElems]),node0_cell)
call discretization_init(microstructureAt,& call discretization_init(materialAt,&
IP_reshaped,& IP_reshaped,&
node0_cell) node0_cell)
@ -172,7 +172,7 @@ end subroutine writeGeometry
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Read mesh from marc input file !> @brief Read mesh from marc input file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt) subroutine inputRead(elem,node0_elem,connectivity_elem,materialAt)
type(tElement), intent(out) :: elem type(tElement), intent(out) :: elem
real(pReal), dimension(:,:), allocatable, intent(out) :: & real(pReal), dimension(:,:), allocatable, intent(out) :: &
@ -180,7 +180,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
integer, dimension(:,:), allocatable, intent(out) :: & integer, dimension(:,:), allocatable, intent(out) :: &
connectivity_elem connectivity_elem
integer, dimension(:), allocatable, intent(out) :: & integer, dimension(:), allocatable, intent(out) :: &
microstructureAt materialAt
integer :: & integer :: &
fileFormatVersion, & fileFormatVersion, &
@ -226,9 +226,9 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,inputFile) connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,inputFile)
call inputRead_microstructure(microstructureAt, & call inputRead_material(materialAt, &
nElems,elem%nNodes,nameElemSet,mapElemSet,& nElems,elem%nNodes,nameElemSet,mapElemSet,&
initialcondTableStyle,inputFile) initialcondTableStyle,inputFile)
end subroutine inputRead end subroutine inputRead
@ -675,13 +675,13 @@ end function inputRead_connectivityElem
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Store microstructure ID !> @brief Store material ID
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine inputRead_microstructure(microstructureAt,& subroutine inputRead_material(materialAt,&
nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileContent) nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileContent)
integer, dimension(:), allocatable, intent(out) :: & integer, dimension(:), allocatable, intent(out) :: &
microstructureAt materialAt
integer, intent(in) :: & integer, intent(in) :: &
nElem, & nElem, &
nNodes, & !< number of nodes per element nNodes, & !< number of nodes per element
@ -696,7 +696,7 @@ subroutine inputRead_microstructure(microstructureAt,&
integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead,l,k,m integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead,l,k,m
allocate(microstructureAt(nElem),source=0) allocate(materialAt(nElem),source=0)
do l = 1, size(fileContent) do l = 1, size(fileContent)
chunkPos = IO_stringPos(fileContent(l)) chunkPos = IO_stringPos(fileContent(l))
@ -715,7 +715,7 @@ subroutine inputRead_microstructure(microstructureAt,&
contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements
do i = 1,contInts(1) do i = 1,contInts(1)
e = mesh_FEM2DAMASK_elem(contInts(1+i)) e = mesh_FEM2DAMASK_elem(contInts(1+i))
microstructureAt(e) = myVal materialAt(e) = myVal
enddo enddo
if (initialcondTableStyle == 0) m = m + 1 if (initialcondTableStyle == 0) m = m + 1
enddo enddo
@ -723,9 +723,9 @@ subroutine inputRead_microstructure(microstructureAt,&
endif endif
enddo enddo
if(any(microstructureAt < 1)) call IO_error(180) if(any(materialAt < 1)) call IO_error(180)
end subroutine inputRead_microstructure end subroutine inputRead_material
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1030,10 +1030,9 @@ pure function IPareaNormal(elem,nElem,connectivity,node)
IPareaNormal(1:3,f,i,e) = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & IPareaNormal(1:3,f,i,e) = math_cross(nodePos(1:3,2) - nodePos(1:3,1), &
nodePos(1:3,3) - nodePos(1:3,1)) nodePos(1:3,3) - nodePos(1:3,1))
case (4) ! 3D 8node case (4) ! 3D 8node
! for this cell type we get the normal of the quadrilateral face as an average of ! Get the normal of the quadrilateral face as the average of four normals of triangular
! four normals of triangular subfaces; since the face consists only of two triangles, ! subfaces. Since the face consists only of two triangles, the sum has to be divided
! the sum has to be divided by two; this whole prcedure tries to compensate for ! by two. This procedure tries to compensate for probable non-planar cell surfaces
! probable non-planar cell surfaces
IPareaNormal(1:3,f,i,e) = 0.0_pReal IPareaNormal(1:3,f,i,e) = 0.0_pReal
do n = 1, m do n = 1, m
IPareaNormal(1:3,f,i,e) = IPareaNormal(1:3,f,i,e) & IPareaNormal(1:3,f,i,e) = IPareaNormal(1:3,f,i,e) &

127
src/material.f90
View File

@ -164,7 +164,7 @@ subroutine material_init(restart)
material_homogenization material_homogenization
character(len=pStringLen) :: sectionName character(len=pStringLen) :: sectionName
print'(/,a)', ' <<<+- material init -+>>>'; flush(6) print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(material_name_phase(phases%length)) allocate(material_name_phase(phases%length))
@ -181,10 +181,10 @@ subroutine material_init(restart)
enddo enddo
call material_parseMicrostructure call material_parseMicrostructure
print*, ' Microstructure parsed' print*, 'Microstructure parsed'
call material_parseHomogenization call material_parseHomogenization
print*, ' Homogenization parsed' print*, 'Homogenization parsed'
if(homogenization_maxNgrains > size(material_phaseAt,1)) call IO_error(148) if(homogenization_maxNgrains > size(material_phaseAt,1)) call IO_error(148)
@ -227,6 +227,7 @@ end subroutine material_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration !> @brief parses the homogenization part from the material configuration
! ToDo: This should be done in homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine material_parseHomogenization subroutine material_parseHomogenization
@ -320,100 +321,78 @@ end subroutine material_parseHomogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine material_parseMicrostructure subroutine material_parseMicrostructure
class(tNode), pointer :: microstructure, & !> pointer to microstructure list class(tNode), pointer :: microstructures, & !> list of microstructures
constituentsInMicrostructure, & !> pointer to a microstructure list item microstructure, & !> microstructure definition
constituents, & !> pointer to constituents list constituents, & !> list of constituents
constituent, & !> pointer to each constituent constituent, & !> constituent definition
phases, & phases, &
homogenization homogenization
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
CounterPhase, & counterPhase, &
CounterHomogenization counterHomogenization
real(pReal), dimension(:,:), allocatable :: &
microstructure_fraction !< vol fraction of each constituent in microstrcuture
real(pReal) :: &
frac
integer :: & integer :: &
e, & e, &
i, & i, &
m, & m, &
c, & c, &
microstructure_maxNconstituents maxNconstituents
real(pReal), dimension(4) :: phase_orientation microstructures => config_material%get('microstructure')
if(any(discretization_microstructureAt > microstructures%length)) &
call IO_error(155,ext_msg='More microstructures requested than found in material.yaml')
homogenization => config_material%get('homogenization') allocate(microstructure_Nconstituents(microstructures%length),source=0)
phases => config_material%get('phase') do m = 1, microstructures%length
microstructure => config_material%get('microstructure') microstructure => microstructures%get(m)
allocate(microstructure_Nconstituents(microstructure%length), source = 0) constituents => microstructure%get('constituents')
if(any(discretization_microstructureAt > microstructure%length)) &
call IO_error(155,ext_msg='More microstructures in geometry than sections in material.yaml')
do m = 1, microstructure%length
constituentsInMicrostructure => microstructure%get(m)
constituents => constituentsInMicrostructure%get('constituents')
microstructure_Nconstituents(m) = constituents%length microstructure_Nconstituents(m) = constituents%length
enddo enddo
maxNconstituents = maxval(microstructure_Nconstituents)
microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
allocate(microstructure_fraction(microstructure_maxNconstituents,microstructure%length), source =0.0_pReal) allocate(material_homogenizationAt(discretization_nElem),source=0)
allocate(material_phaseAt(microstructure_maxNconstituents,discretization_nElem), source =0)
allocate(material_orientation0(microstructure_maxNconstituents,discretization_nIP,discretization_nElem))
allocate(material_homogenizationAt(discretization_nElem))
allocate(material_homogenizationMemberAt(discretization_nIP,discretization_nElem),source=0) allocate(material_homogenizationMemberAt(discretization_nIP,discretization_nElem),source=0)
allocate(material_phaseMemberAt(microstructure_maxNconstituents,discretization_nIP,discretization_nElem),source=0) allocate(material_phaseAt(maxNconstituents,discretization_nElem),source=0)
allocate(material_phaseMemberAt(maxNconstituents,discretization_nIP,discretization_nElem),source=0)
allocate(CounterPhase(phases%length),source=0) allocate(material_orientation0(maxNconstituents,discretization_nIP,discretization_nElem))
allocate(CounterHomogenization(homogenization%length),source=0)
phases => config_material%get('phase')
allocate(counterPhase(phases%length),source=0)
homogenization => config_material%get('homogenization')
allocate(counterHomogenization(homogenization%length),source=0)
do m = 1, microstructure%length do e = 1, discretization_nElem
constituentsInMicrostructure => microstructure%get(m) microstructure => microstructures%get(discretization_microstructureAt(e))
constituents => constituentsInMicrostructure%get('constituents') constituents => microstructure%get('constituents')
material_homogenizationAt(e) = homogenization%getIndex(microstructure%get_asString('homogenization'))
do i = 1, discretization_nIP
counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1
material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e))
enddo
frac = 0.0_pReal
do c = 1, constituents%length do c = 1, constituents%length
constituent => constituents%get(c) constituent => constituents%get(c)
microstructure_fraction(c,m) = constituent%get_asFloat('fraction') frac = frac + constituent%get_asFloat('fraction')
enddo
if (dNeq(sum(microstructure_fraction(:,m)),1.0_pReal)) call IO_error(153,ext_msg='constituent') material_phaseAt(c,e) = phases%getIndex(constituent%get_asString('phase'))
enddo do i = 1, discretization_nIP
counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1
do e = 1, discretization_nElem material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e))
do i = 1, discretization_nIP
constituentsInMicrostructure => microstructure%get(discretization_microstructureAt(e)) call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('orientation',requiredSize=4))
constituents => constituentsInMicrostructure%get('constituents')
do c = 1, constituents%length
constituent => constituents%get(c)
material_phaseAt(c,e) = phases%getIndex(constituent%get_asString('phase'))
phase_orientation = constituent%get_asFloats('orientation')
call material_orientation0(c,i,e)%fromQuaternion(phase_orientation)
enddo enddo
enddo enddo
if (dNeq(frac,1.0_pReal)) call IO_error(153,ext_msg='constituent')
enddo enddo
do e = 1, discretization_nElem
do i = 1, discretization_nIP
constituentsInMicrostructure => microstructure%get(discretization_microstructureAt(e))
material_homogenizationAt(e) = homogenization%getIndex(constituentsInMicrostructure%get_asString('homogenization'))
CounterHomogenization(material_homogenizationAt(e)) = CounterHomogenization(material_homogenizationAt(e)) + 1
material_homogenizationMemberAt(i,e) = CounterHomogenization(material_homogenizationAt(e))
enddo
enddo
do e = 1, discretization_nElem
do i = 1, discretization_nIP
constituentsInMicrostructure => microstructure%get(discretization_microstructureAt(e))
constituents => constituentsInMicrostructure%get('constituents')
do c = 1, constituents%length
CounterPhase(material_phaseAt(c,e)) = &
CounterPhase(material_phaseAt(c,e)) + 1
material_phaseMemberAt(c,i,e) = CounterPhase(material_phaseAt(c,e))
enddo
enddo
enddo
end subroutine material_parseMicrostructure end subroutine material_parseMicrostructure

2
src/math.f90
View File

@ -91,7 +91,7 @@ subroutine math_init
class(tNode), pointer :: & class(tNode), pointer :: &
num_generic num_generic
print'(/,a)', ' <<<+- math init -+>>>'; flush(6) print'(/,a)', ' <<<+- math init -+>>>'; flush(IO_STDOUT)
num_generic => config_numerics%get('generic',defaultVal=emptyDict) num_generic => config_numerics%get('generic',defaultVal=emptyDict)
randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0) randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0)

34
src/mesh/DAMASK_mesh.f90
View File

@ -78,7 +78,7 @@ program DAMASK_mesh
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll call CPFEM_initAll
write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>'; flush(6) print'(/,a)', ' <<<+- DAMASK_mesh init -+>>>'; flush(IO_STDOUT)
!--------------------------------------------------------------------- !---------------------------------------------------------------------
! reading field information from numerics file and do sanity checks ! reading field information from numerics file and do sanity checks
@ -208,30 +208,30 @@ program DAMASK_mesh
errorID = 0 errorID = 0
checkLoadcases: do currentLoadCase = 1, size(loadCases) checkLoadcases: do currentLoadCase = 1, size(loadCases)
write (loadcase_string, '(i0)' ) currentLoadCase write (loadcase_string, '(i0)' ) currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase print'(a,i0)', ' load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
write(6,'(2x,a)') 'drop guessing along trajectory' print'(a)', ' drop guessing along trajectory'
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID) select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label) print'(a)', ' Field '//trim(FIELD_MECH_label)
end select end select
do faceSet = 1, mesh_Nboundaries do faceSet = 1, mesh_Nboundaries
do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) & if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) &
write(6,'(4x,a,i2,a,i2,a,f12.7)') 'Face ', mesh_boundaries(faceSet), & print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), &
' Component ', component, & ' Component ', component, &
' Value ', loadCases(currentLoadCase)%fieldBC(field)% & ' Value ', loadCases(currentLoadCase)%fieldBC(field)% &
componentBC(component)%Value(faceSet) componentBC(component)%Value(faceSet)
enddo enddo
enddo enddo
enddo enddo
write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count
write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs print'(a,i5)', ' increments: ', loadCases(currentLoadCase)%incs
if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', & print'(a,i5)', ' output frequency: ', &
loadCases(currentLoadCase)%outputfrequency loadCases(currentLoadCase)%outputfrequency
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases enddo checkLoadcases
@ -290,8 +290,8 @@ program DAMASK_mesh
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
write(6,'(/,a)') ' ###########################################################################' print'(/,a)', ' ###########################################################################'
write(6,'(1x,a,es12.5,6(a,i0))')& print'(1x,a,es12.5,6(a,i0))',&
'Time', time, & 'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,& '-', stepFraction, '/', subStepFactor**cutBackLevel,&
@ -299,7 +299,7 @@ program DAMASK_mesh
write(incInfo,'(4(a,i0))') & write(incInfo,'(4(a,i0))') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel '-',stepFraction, '/', subStepFactor**cutBackLevel
flush(6) flush(IO_STDOUT)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forward fields ! forward fields
@ -338,7 +338,7 @@ program DAMASK_mesh
cutBack = .False. cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back if (cutBackLevel < maxCutBack) then ! do cut back
write(6,'(/,a)') ' cut back detected' print'(/,a)', ' cut back detected'
cutBack = .True. cutBack = .True.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1 cutBackLevel = cutBackLevel + 1
@ -360,13 +360,13 @@ program DAMASK_mesh
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then if (all(solres(:)%converged)) then
write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' converged' print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged'
else else
write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' NOT converged' print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6) endif; flush(IO_STDOUT)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................' print'(/,a)', ' ... writing results to file ......................................'
call CPFEM_results(totalIncsCounter,time) call CPFEM_results(totalIncsCounter,time)
endif endif
@ -378,7 +378,7 @@ program DAMASK_mesh
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report summary of whole calculation ! report summary of whole calculation
write(6,'(/,a)') ' ###########################################################################' print'(/,a)', ' ###########################################################################'
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
call quit(0) ! no complains ;) call quit(0) ! no complains ;)

2
src/mesh/FEM_quadrature.f90
View File

@ -37,7 +37,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_quadrature_init subroutine FEM_quadrature_init
write(6,'(/,a)') ' <<<+- FEM_quadrature init -+>>>'; flush(6) print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D linear ! 2D linear

8
src/mesh/FEM_utilities.f90
View File

@ -110,7 +110,7 @@ subroutine FEM_utilities_init
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- FEM_utilities init -+>>>' print'(/,a)', ' <<<+- FEM_utilities init -+>>>'
num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2) structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2)
@ -118,11 +118,11 @@ subroutine FEM_utilities_init
debug_mesh => config_debug%get('mesh',defaultVal=emptyList) debug_mesh => config_debug%get('mesh',defaultVal=emptyList)
debugPETSc = debug_mesh%contains('petsc') debugPETSc = debug_mesh%contains('petsc')
if(debugPETSc) write(6,'(3(/,a),/)') & if(debugPETSc) print'(3(/,a),/)', &
' Initializing PETSc with debug options: ', & ' Initializing PETSc with debug options: ', &
trim(PETScDebug), & trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.yaml ' ' add more using the PETSc_Options keyword in numerics.yaml '
flush(6) flush(IO_STDOUT)
call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr) call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
@ -158,7 +158,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' ... evaluating constitutive response ......................................' print'(/,a)', ' ... evaluating constitutive response ......................................'
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field

2
src/mesh/discretization_mesh.f90
View File

@ -83,7 +83,7 @@ subroutine discretization_mesh_init(restart)
num_mesh num_mesh
integer :: integrationOrder !< order of quadrature rule required integer :: integrationOrder !< order of quadrature rule required
write(6,'(/,a)') ' <<<+- discretization_mesh init -+>>>' print'(/,a)', ' <<<+- discretization_mesh init -+>>>'
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
! read numerics parameter ! read numerics parameter

12
src/mesh/mesh_mech_FEM.f90
View File

@ -110,7 +110,7 @@ subroutine FEM_mech_init(fieldBC)
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'; flush(6) print'(/,a)', ' <<<+- FEM_mech init -+>>>'; flush(IO_STDOUT)
!----------------------------------------------------------------------------- !-----------------------------------------------------------------------------
! read numerical parametes and do sanity checks ! read numerical parametes and do sanity checks
@ -318,8 +318,8 @@ type(tSolutionState) function FEM_mech_solution( &
CHKERRQ(ierr) CHKERRQ(ierr)
endif endif
write(6,'(/,a)') ' ===========================================================================' print'(/,a)', ' ==========================================================================='
flush(6) flush(IO_STDOUT)
end function FEM_mech_solution end function FEM_mech_solution
@ -679,12 +679,12 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm
call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr) call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
write(6,'(1/,1x,a,a,i0,a,i0,f0.3)') trim(incInfo), & print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), &
' @ Iteration ',PETScIter,' mechanical residual norm = ', & ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
transpose(P_av)*1.e-6_pReal transpose(P_av)*1.e-6_pReal
flush(6) flush(IO_STDOUT)
end subroutine FEM_mech_converged end subroutine FEM_mech_converged

22
src/parallelization.f90
View File

@ -3,8 +3,8 @@
!> @brief Inquires variables related to parallelization (openMP, MPI) !> @brief Inquires variables related to parallelization (openMP, MPI)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module parallelization module parallelization
use prec use, intrinsic :: ISO_fortran_env, only: &
use, intrinsic :: iso_fortran_env OUTPUT_UNIT
#ifdef PETSc #ifdef PETSc
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
@ -12,6 +12,8 @@ module parallelization
#endif #endif
!$ use OMP_LIB !$ use OMP_LIB
use prec
implicit none implicit none
private private
@ -29,14 +31,14 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine parallelization_init subroutine parallelization_init
integer :: err, typeSize integer :: err, typeSize
!$ integer :: got_env, DAMASK_NUM_THREADS, threadLevel !$ integer :: got_env, DAMASK_NUM_THREADS, threadLevel
!$ character(len=6) NumThreadsString !$ character(len=6) NumThreadsString
#ifdef PETSc #ifdef PETSc
PetscErrorCode :: petsc_err PetscErrorCode :: petsc_err
#else #else
print'(/,a)', ' <<<+- parallelization init -+>>>'; flush(6) print'(/,a)', ' <<<+- parallelization init -+>>>'; flush(OUTPUT_UNIT)
#endif #endif
#ifdef PETSc #ifdef PETSc
@ -69,14 +71,10 @@ subroutine parallelization_init
if (typeSize*8 /= storage_size(0.0_pReal)) error stop 'Mismatch between MPI and DAMASK real' if (typeSize*8 /= storage_size(0.0_pReal)) error stop 'Mismatch between MPI and DAMASK real'
#endif #endif
mainProcess: if (worldrank == 0) then if (worldrank /= 0) then
if (output_unit /= 6) error stop 'STDOUT != 6' close(OUTPUT_UNIT) ! disable output
if (error_unit /= 0) error stop 'STDERR != 0' open(OUTPUT_UNIT,file='/dev/null',status='replace') ! close() alone will leave some temp files in cwd
else mainProcess endif
close(6) ! disable output for non-master processes (open 6 to rank specific file for debug)
open(6,file='/dev/null',status='replace') ! close(6) alone will leave some temp files in cwd
endif mainProcess
!$ call get_environment_variable(name='DAMASK_NUM_THREADS',value=NumThreadsString,STATUS=got_env) !$ call get_environment_variable(name='DAMASK_NUM_THREADS',value=NumThreadsString,STATUS=got_env)
!$ if(got_env /= 0) then !$ if(got_env /= 0) then

2
src/prec.f90
View File

@ -8,7 +8,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module prec module prec
use, intrinsic :: IEEE_arithmetic use, intrinsic :: IEEE_arithmetic
use, intrinsic :: ISO_C_Binding use, intrinsic :: ISO_C_binding
implicit none implicit none
public public

2
src/results.f90
View File

@ -65,7 +65,7 @@ subroutine results_init(restart)
character(len=pStringLen) :: commandLine character(len=pStringLen) :: commandLine
print'(/,a)', ' <<<+- results init -+>>>'; flush(6) print'(/,a)', ' <<<+- results init -+>>>'; flush(IO_STDOUT)
print*, 'Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017' print*, 'Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017'
print*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL print*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL

2
src/rotations.f90
View File

@ -104,7 +104,7 @@ contains
subroutine rotations_init subroutine rotations_init
call quaternions_init call quaternions_init
print'(/,a)', ' <<<+- rotations init -+>>>'; flush(6) print'(/,a)', ' <<<+- rotations init -+>>>'; flush(IO_STDOUT)
print*, 'Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015' print*, 'Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
print*, 'https://doi.org/10.1088/0965-0393/23/8/083501' print*, 'https://doi.org/10.1088/0965-0393/23/8/083501'

2
src/source_damage_anisoBrittle.f90
View File

@ -53,7 +53,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
mySources = source_active('damage_anisoBrittle',source_length) mySources = source_active('damage_anisoBrittle',source_length)
Ninstance = count(mySources) Ninstance = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')

2
src/source_damage_anisoDuctile.f90
View File

@ -47,7 +47,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
mySources = source_active('damage_anisoDuctile',source_length) mySources = source_active('damage_anisoDuctile',source_length)
Ninstance = count(mySources) Ninstance = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')

2
src/source_damage_isoBrittle.f90
View File

@ -43,7 +43,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
mySources = source_active('damage_isoBrittle',source_length) mySources = source_active('damage_isoBrittle',source_length)
Ninstance = count(mySources) Ninstance = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')

2
src/source_damage_isoDuctile.f90
View File

@ -45,7 +45,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
mySources = source_active('damage_isoDuctile',source_length) mySources = source_active('damage_isoDuctile',source_length)
Ninstance = count(mySources) Ninstance = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')

2
src/source_thermal_dissipation.f90
View File

@ -41,7 +41,7 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
mySources = source_active('thermal_dissipation',source_length) mySources = source_active('thermal_dissipation',source_length)
Ninstance = count(mySources) Ninstance = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')

2
src/source_thermal_externalheat.f90
View File

@ -45,7 +45,7 @@ module function source_thermal_externalheat_init(source_length) result(mySources
mySources = source_active('thermal_externalheat',source_length) mySources = source_active('thermal_externalheat',source_length)
Ninstance = count(mySources) Ninstance = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(6) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstance == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')