diff --git a/PRIVATE b/PRIVATE index 86f77da4a..b73dcfe74 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 86f77da4aec6cb52bbcc2724f00a6cf6a7dc6e91 +Subproject commit b73dcfe746eadce92ed0ab08a3bfbb19f5c8eb0a diff --git a/VERSION b/VERSION index a05ebe3cf..a8cc04d22 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha-275-g7801f527f +v3.0.0-alpha-321-g0f6495430 diff --git a/processing/pre/geom_grainGrowth.py b/processing/pre/geom_grainGrowth.py index a573a9fed..b5793f703 100755 --- a/processing/pre/geom_grainGrowth.py +++ b/processing/pre/geom_grainGrowth.py @@ -66,13 +66,13 @@ for name in filenames: grid_original = geom.grid damask.util.croak(geom) - materials = np.tile(geom.materials,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1 - grid = np.array(materials.shape) + material = np.tile(geom.material,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1 + grid = np.array(material.shape) # --- initialize support data --------------------------------------------------------------------- # 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: 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): - interfaceEnergy = np.zeros(materials.shape,dtype=np.float32) + interfaceEnergy = np.zeros(material.shape,dtype=np.float32) for i in (-1,0,1): for j in (-1,0,1): for k in (-1,0,1): # assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood) interfaceEnergy = np.maximum(interfaceEnergy, - getInterfaceEnergy(materials,np.roll(np.roll(np.roll( - materials,i,axis=0), j,axis=1), k,axis=2))) + getInterfaceEnergy(material,np.roll(np.roll(np.roll( + 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 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_indices = True) # want index of closest bulk grain - periodic_materials = np.tile(materials,(3,3,3))[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] # periodically extend the geometry + periodic_material = np.tile(material,(3,3,3))[grid[0]//2:-grid[0]//2, + grid[1]//2:-grid[1]//2, + grid[2]//2:-grid[2]//2] # periodically extend the geometry - materials = periodic_materials[index[0], - index[1], - index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] # extent grains into interface region + material = periodic_material[index[0], + index[1], + index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2, + grid[1]//2:-grid[1]//2, + grid[2]//2:-grid[2]//2] # extent grains into interface region # 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_indices = True) - materials = materials[index[0], + material = material[index[0], index[1], 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 for micro in options.immutable: - immutable += materials_original == micro + immutable += material_original == micro # 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, origin = geom.origin, comments = geom.comments + [scriptID + ' ' + ' '.join(sys.argv[1:])], diff --git a/processing/pre/mentat_spectralBox.py b/processing/pre/mentat_spectralBox.py index 98703219f..7d78cb973 100755 --- a/processing/pre/mentat_spectralBox.py +++ b/processing/pre/mentat_spectralBox.py @@ -100,7 +100,7 @@ def mesh(r,d): #------------------------------------------------------------------------------------------------- -def material(): +def materials(): return [\ "*new_mater standard", "*mater_option general:state:solid", @@ -130,10 +130,10 @@ def geometry(): #------------------------------------------------------------------------------------------------- -def initial_conditions(materials): +def initial_conditions(material): elements = [] element = 0 - for id in materials: + for id in material: element += 1 if len(elements) < id: for i in range(id-len(elements)): @@ -197,14 +197,14 @@ for name in filenames: damask.util.report(scriptName,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 = [\ init(), mesh(geom.grid,geom.size), - material(), + materials(), geometry(), - initial_conditions(materials), + initial_conditions(material), '*identify_sets', '*show_model', '*redraw', diff --git a/processing/pre/seeds_fromDistribution.py b/processing/pre/seeds_fromDistribution.py index 030f3d6fa..48f803d29 100755 --- a/processing/pre/seeds_fromDistribution.py +++ b/processing/pre/seeds_fromDistribution.py @@ -101,8 +101,8 @@ class myThread (threading.Thread): #--- evaluate current seeds file ------------------------------------------------------------------ perturbedGeom = damask.Geom.load_ASCII(perturbedGeomVFile) - myNmaterials = len(np.unique(perturbedGeom.materials)) - currentData=np.bincount(perturbedGeom.materials.ravel())[1:]/points + myNmaterials = len(np.unique(perturbedGeom.material)) + currentData = np.bincount(perturbedGeom.material.ravel())[1:]/points currentError=[] currentHist=[] 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 targetGeom = damask.Geom.load_ASCII(os.path.splitext(os.path.basename(options.target))[0]+'.geom') -nMaterials = len(np.unique(targetGeom.materials)) -targetVolFrac = np.bincount(targetGeom.materials.flatten())/targetGeom.grid.prod().astype(np.float) +nMaterials = len(np.unique(targetGeom.material)) +targetVolFrac = np.bincount(targetGeom.material.flatten())/targetGeom.grid.prod().astype(np.float) target = [] for i in range(1,nMaterials+1): targetHist,targetBins = np.histogram(targetVolFrac,bins=i) #bin boundaries @@ -244,10 +244,10 @@ initialGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+ initialGeomVFile.seek(0) 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') -initialData = np.bincount(initialGeom.materials.flatten())/points +initialData = np.bincount(initialGeom.material.flatten())/points for i in range(nMaterials): initialHist = np.histogram(initialData,bins=target[i]['bins'])[0] 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: - maxSeeds = len(np.unique(initialGeom.materials)) + maxSeeds = len(np.unique(initialGeom.material)) else: maxSeeds = options.maxseeds diff --git a/processing/pre/seeds_fromGeom.py b/processing/pre/seeds_fromGeom.py index de95b439d..b8d74b651 100755 --- a/processing/pre/seeds_fromGeom.py +++ b/processing/pre/seeds_fromGeom.py @@ -47,11 +47,11 @@ for name in filenames: damask.util.report(scriptName,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.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)) 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)\ - .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) diff --git a/processing/pre/seeds_fromPokes.py b/processing/pre/seeds_fromPokes.py index ff9d250fc..887d76392 100755 --- a/processing/pre/seeds_fromPokes.py +++ b/processing/pre/seeds_fromPokes.py @@ -76,7 +76,7 @@ for name in filenames: g[2] = k + offset[2] g %= geom.grid 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.y: g[1] += 1 n += 1 diff --git a/python/damask/_geom.py b/python/damask/_geom.py index a89439028..55569ac14 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -15,13 +15,13 @@ from . import grid_filters class Geom: """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 ---------- - materials : numpy.ndarray + material : numpy.ndarray Material index array (3D). size : list or numpy.ndarray Physical size of the geometry in meter. @@ -31,16 +31,16 @@ class Geom: Comment lines. """ - if len(materials.shape) != 3: - raise ValueError(f'Invalid materials shape {materials.shape}.') - elif materials.dtype not in np.sctypes['float'] + np.sctypes['int']: - raise TypeError(f'Invalid materials data type {materials.dtype}.') + if len(material.shape) != 3: + raise ValueError(f'Invalid material shape {material.shape}.') + elif material.dtype not in np.sctypes['float'] + np.sctypes['int']: + raise TypeError(f'Invalid material data type {material.dtype}.') else: - self.materials = np.copy(materials) + self.material = np.copy(material) - if self.materials.dtype in np.sctypes['float'] and \ - np.all(self.materials == self.materials.astype(int).astype(float)): - self.materials = self.materials.astype(int) + if self.material.dtype in np.sctypes['float'] and \ + np.all(self.material == self.material.astype(int).astype(float)): + self.material = self.material.astype(int) if len(size) != 3 or any(np.array(size) <= 0): raise ValueError(f'Invalid size {size}.') @@ -62,7 +62,7 @@ class Geom: f'size x y z: {util.srepr(self.size, " x ")}', f'origin x y z: {util.srepr(self.origin," ")}', 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.emph( f'# materials: { self.N_materials}')) - if np.nanmax(other.materials) != np.nanmax(self.materials): - message.append(util.delete(f'max material: {np.nanmax(other.materials)}')) - message.append(util.emph( f'max material: {np.nanmax( self.materials)}')) + if np.nanmax(other.material) != np.nanmax(self.material): + message.append(util.delete(f'max material: {np.nanmax(other.material)}')) + message.append(util.emph( f'max material: {np.nanmax( self.material)}')) return util.return_message(message) @property def grid(self): - return np.asarray(self.materials.shape) + return np.asarray(self.material.shape) @property def N_materials(self): - return np.unique(self.materials).size + return np.unique(self.material).size @staticmethod @@ -160,7 +160,7 @@ class Geom: else: comments.append(line.strip()) - materials = np.empty(grid.prod()) # initialize as flat array + material = np.empty(grid.prod()) # initialize as flat array i = 0 for line in content[header_length:]: items = line.split('#')[0].split() @@ -172,16 +172,16 @@ class Geom: abs(int(items[2])-int(items[0]))+1,dtype=float) 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) if i != grid.prod(): raise TypeError(f'Invalid file: expected {grid.prod()} entries, found {i}') - if not np.any(np.mod(materials,1) != 0.0): # no float present - materials = materials.astype('int') + if not np.any(np.mod(material,1) != 0.0): # no float present + 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 @@ -200,9 +200,11 @@ class Geom: comments = v.get_comments() grid = np.array(v.vtk_data.GetDimensions())-1 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 @@ -210,7 +212,7 @@ class Geom: return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights) @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. @@ -224,7 +226,7 @@ class Geom: Position of the seed points in meter. All points need to lay within the box. weights : numpy.ndarray of shape (seeds.shape[0]) 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 consecutively numbered. 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]) pool.close() pool.join() - materials_ = np.array(result.get()) + material_ = np.array(result.get()) if periodic: - materials_ = materials_.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_.reshape(grid*3) + material_ = material_[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] else: - materials_ = materials_.reshape(grid) + material_ = material_.reshape(grid) - geom = Geom(materials = materials_+1, - size = size, - comments = util.execution_stamp('Geom','from_Laguerre_tessellation'), + geom = Geom(material = material_+1, + size = size, + comments = util.execution_stamp('Geom','from_Laguerre_tessellation'), ) - if materials is not None: - geom = geom.substitute(np.arange(seeds.shape[0])+1,materials) + if material is not None: + geom = geom.substitute(np.arange(seeds.shape[0])+1,material) geom.comments = geom.comments[:-1] return geom @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. @@ -278,7 +280,7 @@ class Geom: Physical size of the geometry in meter. seeds : numpy.ndarray of shape (:,3) 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 consecutively numbered. periodic : Boolean, optional @@ -287,14 +289,14 @@ class Geom: """ coords = grid_filters.cell_coord0(grid,size).reshape(-1,3) 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, - size = size, - comments = util.execution_stamp('Geom','from_Voronoi_tessellation'), + geom = Geom(material = material_.reshape(grid)+1, + size = size, + comments = util.execution_stamp('Geom','from_Voronoi_tessellation'), ) - if materials is not None: - geom = geom.substitute(np.arange(seeds.shape[0])+1,materials) + if material is not None: + geom = geom.substitute(np.arange(seeds.shape[0])+1,material) geom.comments = geom.comments[:-1] return geom @@ -327,10 +329,10 @@ class Geom: plain = not compress if plain: - format_string = '%g' if self.materials.dtype in np.sctypes['float'] else \ - '%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.materials))))) + format_string = '%g' if self.material.dtype in np.sctypes['float'] else \ + '%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.material))))) 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='') else: try: @@ -341,7 +343,7 @@ class Geom: compressType = None former = start = -1 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): compressType = 'to' reps += 1 @@ -386,7 +388,7 @@ class Geom: """ 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.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) 1 gives a sphere (|x|^(2^1) + |y|^(2^1) + |z|^(2^1) < 1) 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 Rotation of primitive. Defaults to no rotation. inverse : Boolean, optional @@ -444,12 +446,12 @@ class Geom: if periodic: # translate back to center 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_), - size = self.size, - origin = self.origin, - comments = self.comments+[util.execution_stamp('Geom','add_primitive')], + return Geom(material = np.where(np.logical_not(mask) if inverse else mask, self.material,fill_), + size = self.size, + origin = self.origin, + 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.') limits = [None,None] if reflect else [-2,0] - ms = self.materials.copy() + mat = self.material.copy() 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: - 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: - 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, - size = self.size/self.grid*np.asarray(ms.shape), - origin = self.origin, - comments = self.comments+[util.execution_stamp('Geom','mirror')], + return Geom(material = mat, + size = self.size/self.grid*np.asarray(mat.shape), + origin = self.origin, + comments = self.comments+[util.execution_stamp('Geom','mirror')], ) @@ -502,12 +504,12 @@ class Geom: if not set(directions).issubset(valid): 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, - size = self.size, - origin = self.origin, - comments = self.comments+[util.execution_stamp('Geom','flip')], + return Geom(material = mat, + size = self.size, + origin = self.origin, + comments = self.comments+[util.execution_stamp('Geom','flip')], ) @@ -523,17 +525,17 @@ class Geom: Assume geometry to be periodic. Defaults to True. """ - return Geom(materials = ndimage.interpolation.zoom( - self.materials, + return Geom(material = ndimage.interpolation.zoom( + self.material, grid/self.grid, - output=self.materials.dtype, + output=self.material.dtype, order=0, mode=('wrap' if periodic else 'nearest'), prefilter=False ), - size = self.size, - origin = self.origin, - comments = self.comments+[util.execution_stamp('Geom','scale')], + size = self.size, + origin = self.origin, + comments = self.comments+[util.execution_stamp('Geom','scale')], ) @@ -559,29 +561,29 @@ class Geom: else: return me - return Geom(materials = ndimage.filters.generic_filter( - self.materials, + return Geom(material = ndimage.filters.generic_filter( + self.material, mostFrequent, size=(stencil if selection is None else stencil//2*2+1,)*3, mode=('wrap' if periodic else 'nearest'), extra_keywords=dict(selection=selection), - ).astype(self.materials.dtype), - size = self.size, - origin = self.origin, - comments = self.comments+[util.execution_stamp('Geom','clean')], + ).astype(self.material.dtype), + size = self.size, + origin = self.origin, + comments = self.comments+[util.execution_stamp('Geom','clean')], ) def renumber(self): """Renumber sorted material indices to 1,...,N.""" - renumbered = np.empty(self.grid,dtype=self.materials.dtype) - for i, oldID in enumerate(np.unique(self.materials)): - renumbered = np.where(self.materials == oldID, i+1, renumbered) + renumbered = np.empty(self.grid,dtype=self.material.dtype) + for i, oldID in enumerate(np.unique(self.material)): + renumbered = np.where(self.material == oldID, i+1, renumbered) - return Geom(materials = renumbered, - size = self.size, - origin = self.origin, - comments = self.comments+[util.execution_stamp('Geom','renumber')], + return Geom(material = renumbered, + size = self.size, + origin = self.origin, + comments = self.comments+[util.execution_stamp('Geom','renumber')], ) @@ -594,32 +596,32 @@ class Geom: R : damask.Rotation Rotation to apply to the geometry. 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 - dtype = float if np.isnan(fill) or int(fill) != fill or self.materials.dtype==np.float else int + if fill is None: fill = np.nanmax(self.material) + 1 + dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int 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'') # 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)]): - 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) - 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° - 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: - 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, - size = self.size/self.grid*np.asarray(materials_in.shape), - origin = origin, - comments = self.comments+[util.execution_stamp('Geom','rotate')], + return Geom(material = material_in, + size = self.size/self.grid*np.asarray(material_in.shape), + origin = origin, + comments = self.comments+[util.execution_stamp('Geom','rotate')], ) @@ -634,12 +636,12 @@ class Geom: offset : numpy.ndarray of shape (3) Offset (measured in grid points) from old to new geometry [0,0,0]. 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 fill is None: fill = np.nanmax(self.materials) + 1 - dtype = float if int(fill) != fill or self.materials.dtype in np.sctypes['float'] else int + if fill is None: fill = np.nanmax(self.material) + 1 + 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) @@ -648,35 +650,35 @@ class Geom: ll = np.clip(-offset, 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, - size = self.size/self.grid*np.asarray(canvas.shape), - origin = self.origin+offset*self.size/self.grid, - comments = self.comments+[util.execution_stamp('Geom','canvas')], + return Geom(material = canvas, + size = self.size/self.grid*np.asarray(canvas.shape), + origin = self.origin+offset*self.size/self.grid, + 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. Parameters ---------- - from_materials : iterable of ints + from_material : iterable of ints Material indices to be substituted. - to_materials : iterable of ints + to_material : iterable of ints New material indices. """ - substituted = self.materials.copy() - for from_ms,to_ms in zip(from_materials,to_materials): - substituted[self.materials==from_ms] = to_ms + substituted = self.material.copy() + for from_ms,to_ms in zip(from_material,to_material): + substituted[self.material==from_ms] = to_ms - return Geom(materials = substituted, - size = self.size, - origin = self.origin, - comments = self.comments+[util.execution_stamp('Geom','substitute')], + return Geom(material = substituted, + size = self.size, + origin = self.origin, + comments = self.comments+[util.execution_stamp('Geom','substitute')], ) @@ -695,7 +697,7 @@ class Geom: Defaults to 1. offset : int, optional Offset (positive or negative) to tag material indices, - defaults to materials.max() + 1. + defaults to material.max() + 1. trigger : list of ints, optional List of material indices that trigger a change. Defaults to [], meaning that any different neighbor triggers a change. @@ -714,15 +716,15 @@ class Geom: trigger = list(trigger) return np.any(np.in1d(stencil,np.array(trigger))) - offset_ = np.nanmax(self.materials) if offset is None else offset - mask = ndimage.filters.generic_filter(self.materials, + offset_ = np.nanmax(self.material) if offset is None else offset + mask = ndimage.filters.generic_filter(self.material, tainted_neighborhood, size=1+2*vicinity, mode='wrap' if periodic else 'nearest', extra_keywords={'trigger':trigger}) - return Geom(materials = np.where(mask, self.materials + offset_,self.materials), - size = self.size, - origin = self.origin, - comments = self.comments+[util.execution_stamp('Geom','vicinity_offset')], + return Geom(material = np.where(mask, self.material + offset_,self.material), + size = self.size, + origin = self.origin, + comments = self.comments+[util.execution_stamp('Geom','vicinity_offset')], ) diff --git a/python/damask/util.py b/python/damask/util.py index 6432ed4e1..43e86314e 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -117,6 +117,7 @@ def execute(cmd, initialPath = os.getcwd() myEnv = os.environ if env is None else env os.chdir(wd) + print(f"executing '{cmd}' in '{wd}'") process = subprocess.Popen(shlex.split(cmd), stdout = subprocess.PIPE, stderr = subprocess.PIPE, @@ -128,7 +129,7 @@ def execute(cmd, stdout = stdout.decode('utf-8').replace('\x08','') stderr = stderr.decode('utf-8').replace('\x08','') 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 @@ -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 = 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))]): raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?') diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 957bc69e6..a91365727 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -11,8 +11,8 @@ from damask import util def geom_equal(a,b): - return np.all(a.materials == b.materials) and \ - np.all(a.grid == b.grid) and \ + return np.all(a.material == b.material) and \ + np.all(a.grid == b.grid) and \ np.allclose(a.size, b.size) and \ str(a.diff(b)) == str(b.diff(a)) @@ -38,7 +38,7 @@ class TestGeom: 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)) != '' @@ -93,28 +93,28 @@ class TestGeom: def test_invalid_size(self,default): with pytest.raises(ValueError): - Geom(default.materials[1:,1:,1:], + Geom(default.material[1:,1:,1:], size=np.ones(2)) def test_invalid_origin(self,default): with pytest.raises(ValueError): - Geom(default.materials[1:,1:,1:], + Geom(default.material[1:,1:,1:], size=np.ones(3), origin=np.ones(4)) def test_invalid_materials_shape(self,default): - materials = np.ones((3,3)) + material = np.ones((3,3)) with pytest.raises(ValueError): - Geom(materials, + Geom(material, size=np.ones(3)) 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): - Geom(materials) + Geom(material) @pytest.mark.parametrize('directions,reflect',[ @@ -205,10 +205,10 @@ class TestGeom: def test_renumber(self,default): - materials = default.materials.copy() - for m in np.unique(materials): - materials[materials==m] = materials.max() + np.random.randint(1,30) - modified = Geom(materials, + material = default.material.copy() + for m in np.unique(material): + material[material==m] = material.max() + np.random.randint(1,30) + modified = Geom(material, default.size, default.origin) assert not geom_equal(modified,default) @@ -218,13 +218,13 @@ class TestGeom: def test_substitute(self,default): offset = np.random.randint(1,500) - modified = Geom(default.materials + offset, + modified = Geom(default.material + offset, default.size, default.origin) assert not geom_equal(modified,default) assert geom_equal(default, - modified.substitute(np.arange(default.materials.max())+1+offset, - np.arange(default.materials.max())+1)) + modified.substitute(np.arange(default.material.max())+1+offset, + 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]), @@ -251,7 +251,7 @@ class TestGeom: grid = default.grid grid_add = np.random.randint(0,30,(3)) 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), @@ -271,7 +271,7 @@ class TestGeom: s = np.random.random(3)+.5 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) - 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)), @@ -308,14 +308,14 @@ class TestGeom: 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]) def test_vicinity_offset_invariant(self,default,periodic): - offset = default.vicinity_offset(trigger=[default.materials.max()+1, - default.materials.min()-1]) - assert np.all(offset.materials==default.materials) + offset = default.vicinity_offset(trigger=[default.material.max()+1, + default.material.min()-1]) + assert np.all(offset.material==default.material) @pytest.mark.parametrize('periodic',[True,False]) @@ -338,7 +338,7 @@ class TestGeom: ms = np.random.randint(1, N_seeds+1) weights[ms-1] = np.random.random() 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']) @@ -346,10 +346,10 @@ class TestGeom: grid = np.random.randint(5,10,3)*2 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]))) - materials = np.ones(grid) - materials[:,grid[1]//2:,:] = 2 + material = np.ones(grid) + material[:,grid[1]//2:,:] = 2 if approach == 'Laguerre': geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),periodic=np.random.random()>0.5) elif approach == 'Voronoi': 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) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index d0fc6413e..05255e2a1 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -106,7 +106,7 @@ subroutine CPFEM_init num_commercialFEM, & 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_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_dcsdE: ', shape(CPFEM_dcsdE) print'(a32,1x,6(i8,1x),/)', 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood) - flush(6) + flush(IO_STDOUT) endif 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 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 - flush(6) + flush(IO_STDOUT) endif endif diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index c87f361c2..fed43ba78 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -76,7 +76,7 @@ end subroutine CPFEM_initAll !-------------------------------------------------------------------------------------------------- subroutine CPFEM_init - print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(6) + print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT) if (interface_restartInc > 0) call crystallite_restartRead diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index 46ecd4da0..d34034d0d 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -14,7 +14,7 @@ #define PETSC_MINOR_MAX 13 module DAMASK_interface - use, intrinsic :: iso_fortran_env + use, intrinsic :: ISO_fortran_env use PETScSys @@ -82,7 +82,7 @@ subroutine 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 #ifdef DEBUG @@ -101,8 +101,8 @@ subroutine DAMASK_interface_init #endif print*, achar(27)//'[0m' - print'(a)', ' Roters et al., Computational Materials Science 158:420–478, 2019' - print'(a)', ' https://doi.org/10.1016/j.commatsci.2018.04.030' + print*, 'Roters et al., Computational Materials Science 158:420–478, 2019' + print*, 'https://doi.org/10.1016/j.commatsci.2018.04.030' print'(/,a)', ' Version: '//DAMASKVERSION @@ -373,7 +373,7 @@ function makeRelativePath(a,b) a_cleaned = rectifyPath(trim(a)//'/') 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) == '/') posLastCommonSlash = i enddo @@ -395,7 +395,7 @@ subroutine catchSIGTERM(signal) bind(C) integer(C_INT), value :: signal 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 @@ -420,7 +420,7 @@ subroutine catchSIGUSR1(signal) bind(C) integer(C_INT), value :: signal 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 @@ -445,7 +445,7 @@ subroutine catchSIGUSR2(signal) bind(C) integer(C_INT), value :: signal 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 diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index b6f9436ff..ea7430c6b 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -30,7 +30,7 @@ module DAMASK_interface use prec #if __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & + use, intrinsic :: ISO_fortran_env, only: & compiler_version, & compiler_options #endif diff --git a/src/IO.f90 b/src/IO.f90 index 6da7a10c8..f9e708ead 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -6,6 +6,10 @@ !> @brief input/output functions !-------------------------------------------------------------------------------------------------- module IO + use, intrinsic :: ISO_fortran_env, only: & + IO_STDOUT => OUTPUT_UNIT, & + IO_STDERR => ERROR_UNIT + use prec implicit none @@ -16,7 +20,7 @@ module IO character, parameter, public :: & IO_EOL = new_line('DAMASK'), & !< end of line character IO_COMMENT = '#' - character(len=*), parameter, private :: & + character(len=*), parameter :: & IO_DIVIDER = '───────────────────'//& '───────────────────'//& '───────────────────'//& @@ -37,7 +41,8 @@ module IO IO_stringAsFloat, & IO_stringAsBool, & IO_error, & - IO_warning + IO_warning, & + IO_STDOUT contains @@ -47,7 +52,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine IO_init - print'(/,a)', ' <<<+- IO init -+>>>'; flush(6) + print'(/,a)', ' <<<+- IO init -+>>>'; flush(IO_STDOUT) call selfTest @@ -538,29 +543,29 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) end select !$OMP CRITICAL (write2out) - write(0,'(/,a)') ' ┌'//IO_DIVIDER//'┐' - write(0,'(a,24x,a,40x,a)') ' │','error', '│' - write(0,'(a,24x,i3,42x,a)') ' │',error_ID, '│' - write(0,'(a)') ' ├'//IO_DIVIDER//'┤' + write(IO_STDERR,'(/,a)') ' ┌'//IO_DIVIDER//'┐' + write(IO_STDERR,'(a,24x,a,40x,a)') ' │','error', '│' + write(IO_STDERR,'(a,24x,i3,42x,a)') ' │',error_ID, '│' + write(IO_STDERR,'(a)') ' ├'//IO_DIVIDER//'┤' 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)' - write(0,formatString) '│ ',trim(msg), '│' + write(IO_STDERR,formatString) '│ ',trim(msg), '│' if (present(ext_msg)) then 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)' - write(0,formatString) '│ ',trim(ext_msg), '│' + write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│' endif 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)) & - write(0,'(a19,1x,i9,44x,a3)') ' │ at IP ',ip, '│' + write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at IP ',ip, '│' 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)) & - write(0,'(a19,1x,i9,44x,a3)') ' │ at instance ',instance, '│' - write(0,'(a,69x,a)') ' │', '│' - write(0,'(a)') ' └'//IO_DIVIDER//'┘' - flush(0) + write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at instance ',instance, '│' + write(IO_STDERR,'(a,69x,a)') ' │', '│' + write(IO_STDERR,'(a)') ' └'//IO_DIVIDER//'┘' + flush(IO_STDERR) call quit(9000+error_ID) !$OMP END CRITICAL (write2out) @@ -623,27 +628,27 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) end select !$OMP CRITICAL (write2out) - write(6,'(/,a)') ' ┌'//IO_DIVIDER//'┐' - write(6,'(a,24x,a,38x,a)') ' │','warning', '│' - write(6,'(a,24x,i3,42x,a)') ' │',warning_ID, '│' - write(6,'(a)') ' ├'//IO_DIVIDER//'┤' + write(IO_STDERR,'(/,a)') ' ┌'//IO_DIVIDER//'┐' + write(IO_STDERR,'(a,24x,a,38x,a)') ' │','warning', '│' + write(IO_STDERR,'(a,24x,i3,42x,a)') ' │',warning_ID, '│' + write(IO_STDERR,'(a)') ' ├'//IO_DIVIDER//'┤' 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)' - write(6,formatString) '│ ',trim(msg), '│' + write(IO_STDERR,formatString) '│ ',trim(msg), '│' if (present(ext_msg)) then 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)' - write(6,formatString) '│ ',trim(ext_msg), '│' + write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│' endif 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)) & - write(6,'(a19,1x,i9,44x,a3)') ' │ at IP ',ip, '│' + write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at IP ',ip, '│' if (present(g)) & - write(6,'(a19,1x,i9,44x,a3)') ' │ at constituent',g, '│' - write(6,'(a,69x,a)') ' │', '│' - write(6,'(a)') ' └'//IO_DIVIDER//'┘' - flush(6) + write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at constituent',g, '│' + write(IO_STDERR,'(a,69x,a)') ' │', '│' + write(IO_STDERR,'(a)') ' └'//IO_DIVIDER//'┘' + flush(IO_STDERR) !$OMP END CRITICAL (write2out) end subroutine IO_warning diff --git a/src/base64.f90 b/src/base64.f90 index 3a59b7049..2f91334b7 100644 --- a/src/base64.f90 +++ b/src/base64.f90 @@ -27,7 +27,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine base64_init - print'(/,a)', ' <<<+- base64 init -+>>>'; flush(6) + print'(/,a)', ' <<<+- base64 init -+>>>'; flush(IO_STDOUT) call selfTest diff --git a/src/config.f90 b/src/config.f90 index 50d4e96e8..b10edf013 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -35,7 +35,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine config_init - print'(/,a)', ' <<<+- config init -+>>>'; flush(6) + print'(/,a)', ' <<<+- config init -+>>>'; flush(IO_STDOUT) call parse_material call parse_numerics @@ -59,7 +59,7 @@ subroutine parse_material inquire(file=fname,exist=fileExists) if(.not. fileExists) call IO_error(100,ext_msg=fname) endif - print*, 'reading '//fname; flush(6) + print*, 'reading '//fname; flush(IO_STDOUT) config_material => YAML_parse_file(fname) end subroutine parse_material @@ -75,7 +75,7 @@ subroutine parse_numerics config_numerics => emptyDict inquire(file='numerics.yaml', exist=fexist) if (fexist) then - print*, 'reading numerics.yaml'; flush(6) + print*, 'reading numerics.yaml'; flush(IO_STDOUT) config_numerics => YAML_parse_file('numerics.yaml') endif @@ -92,7 +92,7 @@ subroutine parse_debug config_debug => emptyDict inquire(file='debug.yaml', exist=fexist) fileExists: if (fexist) then - print*, 'reading debug.yaml'; flush(6) + print*, 'reading debug.yaml'; flush(IO_STDOUT) config_debug => YAML_parse_file('debug.yaml') endif fileExists diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 3b342e3e6..62846359d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -446,7 +446,7 @@ subroutine constitutive_init call damage_init call thermal_init - print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(6) + print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) constitutive_source_maxSizeDotState = 0 PhaseLoop2:do p = 1,phases%length diff --git a/src/constitutive_plastic_disloTungsten.f90 b/src/constitutive_plastic_disloTungsten.f90 index 2bf4fd48e..98a9ce2f6 100644 --- a/src/constitutive_plastic_disloTungsten.f90 +++ b/src/constitutive_plastic_disloTungsten.f90 @@ -100,7 +100,7 @@ module function plastic_disloTungsten_init() result(myPlasticity) myPlasticity = plastic_active('disloTungsten') Ninstance = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016' diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index a25815899..4890316b8 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -147,7 +147,7 @@ module function plastic_dislotwin_init() result(myPlasticity) myPlasticity = plastic_active('dislotwin') Ninstance = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index c78d497d6..477938145 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -71,7 +71,7 @@ module function plastic_isotropic_init() result(myPlasticity) myPlasticity = plastic_active('isotropic') Ninstance = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index fe17b090e..55e482db6 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -83,7 +83,7 @@ module function plastic_kinehardening_init() result(myPlasticity) myPlasticity = plastic_active('kinehardening') Ninstance = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return allocate(param(Ninstance)) diff --git a/src/constitutive_plastic_none.f90 b/src/constitutive_plastic_none.f90 index d62a798cc..ab5f32d3f 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -35,7 +35,7 @@ module function plastic_none_init() result(myPlasticity) enddo Ninstance = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return do p = 1, phases%length diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index c96301691..21523161b 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -189,7 +189,7 @@ module function plastic_nonlocal_init() result(myPlasticity) myPlasticity = plastic_active('nonlocal') Ninstance = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) then call geometry_plastic_nonlocal_disable 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*, '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(state(Ninstance)) @@ -741,10 +741,10 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) if (debugConstitutive%extensive & .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)& .or. .not. debugConstitutive%selective)) then - write(6,'(/,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) - write(6,'(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,i8,1x,i2,1x,i1,/)', '<< CONST >> nonlocal_microstructure at el ip ',el,ip + print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,of) + print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,of)*1e-6 + print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,of)*1e-6 endif #endif @@ -958,8 +958,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el) if (debugConstitutive%extensive & .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)& .or. .not. debugConstitutive%selective)) then - write(6,'(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,/,8(12x,12(e12.5,1x),/))', '<< CONST >> dislocation remobilization', deltaRhoRemobilization(:,1:8) + print'(a,/,10(12x,12(e12.5,1x),/),/)', '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress endif #endif @@ -1047,8 +1047,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & if (debugConstitutive%basic & .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip) & .or. .not. debugConstitutive%selective)) then - write(6,'(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,/,10(12x,12(e12.5,1x),/))', '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip + print'(a,/,4(12x,12(e12.5,1x),/))', '<< CONST >> gdot / 1/s',gdot 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 #ifdef DEBUG if (debugConstitutive%extensive) then - write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip - write(6,'(a)') '<< CONST >> enforcing cutback !!!' + print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip + print'(a)', '<< CONST >> enforcing cutback !!!' endif #endif 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) #ifdef DEBUG if (debugConstitutive%extensive) then - write(6,'(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,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip + print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & maxval(abs(v0), abs(gdot) > 0.0_pReal & .and. prm%CFLfactor * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index b30a4d9df..8cf63a54d 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -92,7 +92,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) myPlasticity = plastic_active('phenopowerlaw') Ninstance = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return allocate(param(Ninstance)) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 2ee85024b..392dba277 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -294,7 +294,7 @@ subroutine crystallite_init print'(a42,1x,i10)', ' # of elements: ', eMax print'(a42,1x,i10)', ' # of integration points/element: ', iMax print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax - flush(6) + flush(IO_STDOUT) endif #endif @@ -1561,7 +1561,7 @@ subroutine crystallite_restartWrite integer(HID_T) :: fileHandle, groupHandle 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' fileHandle = HDF5_openFile(fileName,'a') diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 3588010b2..af2532184 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -49,7 +49,7 @@ subroutine damage_local_init homog, & homogDamage - print'(/,a)', ' <<<+- damage_local init -+>>>'; flush(6) + print'(/,a)', ' <<<+- damage_local init -+>>>'; flush(IO_STDOUT) !---------------------------------------------------------------------------------------------- ! read numerics parameter and do sanity check diff --git a/src/element.f90 b/src/element.f90 index 3ccc2fb78..722a7fd96 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -922,7 +922,7 @@ subroutine tElement_init(self,elemType) 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*, ' geom type: ',self%geomType diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index a5001de44..8158996d6 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -99,10 +99,10 @@ program DAMASK_grid ! init DAMASK (all modules) 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' - write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' + print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' + print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' !-------------------------------------------------------------------------------------------------- ! initialize field solver information @@ -263,56 +263,56 @@ program DAMASK_grid reportAndCheck: if (worldrank == 0) then write (loadcase_string, '(i0)' ) currentLoadCase - write(6,'(/,1x,a,i0)') 'load case: ', currentLoadCase + print'(/,a,i0)', ' load case: ', currentLoadCase if (.not. newLoadCase%followFormerTrajectory) & - write(6,'(2x,a)') 'drop guessing along trajectory' + print*, ' drop guessing along trajectory' if (newLoadCase%deformation%myType == 'l') then do j = 1, 3 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 enddo - write(6,'(2x,a)') 'velocity gradient:' + print*, ' velocity gradient:' 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 - write(6,'(2x,a)') 'deformation gradient rate:' + print*, ' deformation gradient rate:' endif do i = 1, 3; do j = 1, 3 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 - write(6,'(2x,12a)',advance='no') ' * ' + write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' endif - enddo; write(6,'(/)',advance='no') + enddo; write(IO_STDOUT,'(/)',advance='no') enddo if (any(newLoadCase%stress%maskLogical .eqv. & newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) & .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 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 - write(6,'(2x,12a)',advance='no') ' * ' + write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' endif - enddo; write(6,'(/)',advance='no') + enddo; write(IO_STDOUT,'(/)',advance='no') enddo if (any(abs(matmul(newLoadCase%rot%asMatrix(), & transpose(newLoadCase%rot%asMatrix()))-math_I3) > & reshape(spread(tol_math_check,1,9),[ 3,3]))) errorID = 846 ! given rotation matrix contains strain 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()) 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 - write(6,'(2x,a,i0)') 'increments: ', newLoadCase%incs + print'(a,i0)', ' increments: ', newLoadCase%incs 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 < 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 endif reportAndCheck loadCases = [loadCases,newLoadCase] ! load case is ok, append it @@ -341,9 +341,8 @@ program DAMASK_grid writeHeader: if (interface_restartInc < 1) then open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file - if (debug_grid%contains('basic')) & - write(6,'(/,a)') ' header of statistics file written out' - flush(6) + if (debug_grid%contains('basic')) print'(/,a)', ' header of statistics file written out' + flush(IO_STDOUT) else writeHeader open(newunit=statUnit,file=trim(getSolverJobName())//& '.sta',form='FORMATTED', position='APPEND', status='OLD') @@ -351,7 +350,7 @@ program DAMASK_grid endif 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) endif writeUndeformed @@ -397,8 +396,8 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! report begin of new step - write(6,'(/,a)') ' ###########################################################################' - write(6,'(1x,a,es12.5,6(a,i0))') & + print'(/,a)', ' ###########################################################################' + print'(1x,a,es12.5,6(a,i0))', & 'Time', time, & 's: Increment ', inc,'/',loadCases(currentLoadCase)%incs,& '-', stepFraction,'/',subStepFactor**cutBackLevel,& @@ -406,7 +405,7 @@ program DAMASK_grid write(incInfo,'(4(a,i0))') & 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& '-', stepFraction,'/',subStepFactor**cutBackLevel - flush(6) + flush(IO_STDOUT) !-------------------------------------------------------------------------------------------------- ! forward fields @@ -475,7 +474,7 @@ program DAMASK_grid cutBackLevel = cutBackLevel + 1 time = time - timeinc ! rewind time timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep - write(6,'(/,a)') ' cutting back ' + print'(/,a)', ' cutting back ' else ! no more options to continue call IO_warning(850) if (worldrank == 0) close(statUnit) @@ -487,14 +486,14 @@ program DAMASK_grid cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc if (all(solres(:)%converged)) then - write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' converged' + print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged' else - write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' NOT converged' - endif; flush(6) + print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged' + endif; flush(IO_STDOUT) if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency - write(6,'(1/,a)') ' ... writing results to file ......................................' - flush(6) + print'(1/,a)', ' ... writing results to file ......................................' + flush(IO_STDOUT) call CPFEM_results(totalIncsCounter,time) endif if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then @@ -510,7 +509,7 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation - write(6,'(/,a)') ' ###########################################################################' + print'(/,a)', ' ###########################################################################' if (worldrank == 0) close(statUnit) call quit(0) ! no complains ;) diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 276dd566f..27df0acd5 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -56,7 +56,7 @@ subroutine discretization_grid_init(restart) myGrid !< domain grid of this process integer, dimension(:), allocatable :: & - microstructureAt + materialAt integer :: & j, & @@ -65,12 +65,12 @@ subroutine discretization_grid_init(restart) integer(C_INTPTR_T) :: & 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 - call readVTR(grid,geomSize,origin,microstructureAt) + call readVTR(grid,geomSize,origin,materialAt) else - call readGeom(grid,geomSize,origin,microstructureAt) + call readGeom(grid,geomSize,origin,materialAt) endif print'(/,a,3(i12 ))', ' grid a b c: ', grid @@ -102,15 +102,14 @@ subroutine discretization_grid_init(restart) !-------------------------------------------------------------------------------------------------- ! general discretization - microstructureAt = microstructureAt(product(grid(1:2))*grid3Offset+1: & - product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI + materialAt = materialAt(product(grid(1:2))*grid3Offset+1:product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI - call discretization_init(microstructureAt, & + call discretization_init(materialAt, & IPcoordinates0(myGrid,mySize,grid3Offset), & Nodes0(myGrid,mySize,grid3Offset),& - merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write bottom layer - (grid(1)+1) * (grid(2)+1) * grid3,& ! do not write bottom layer (is top of rank-1) - worldrank<1)) + merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write top layer... + (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process + worldrank+1==worldsize)) 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 @@ -147,7 +146,7 @@ end subroutine discretization_grid_init !> @details important variables have an implicit "save" attribute. Therefore, this function is ! supposed to be called only once! !-------------------------------------------------------------------------------------------------- -subroutine readGeom(grid,geomSize,origin,microstructure) +subroutine readGeom(grid,geomSize,origin,material) integer, dimension(3), intent(out) :: & grid ! grid (across all processes!) @@ -155,7 +154,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure) geomSize, & ! size (across all processes!) origin ! origin (across all processes!) integer, dimension(:), intent(out), allocatable :: & - microstructure + material character(len=:), allocatable :: rawData character(len=65536) :: line @@ -167,7 +166,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure) startPos, endPos, & myStat, & l, & !< line counter - c, & !< counter for # microstructures in line + c, & !< counter for # materials in line o, & !< order of "to" packing e, & !< "element", i.e. spectral collocation point i, j @@ -266,7 +265,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure) if(any(geomSize < 0.0_pReal)) & 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 @@ -281,18 +280,18 @@ subroutine readGeom(grid,geomSize,origin,microstructure) noCompression: if (chunkPos(1) /= 3) then 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 compression: if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'of') then 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 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)) - 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 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 noCompression @@ -308,7 +307,7 @@ end subroutine readGeom !> @brief Parse vtk rectilinear grid (.vtr) !> @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) :: & grid ! grid (across all processes!) @@ -316,7 +315,7 @@ subroutine readVTR(grid,geomSize,origin,microstructure) geomSize, & ! size (across all processes!) origin ! origin (across all processes!) integer, dimension(:), intent(out), allocatable :: & - microstructure + material character(len=:), allocatable :: fileContent, dataType, headerType logical :: inFile,inGrid,gotCoordinates,gotCellData,compressed @@ -364,11 +363,9 @@ subroutine readVTR(grid,geomSize,origin,microstructure) else if(index(fileContent(startPos:endPos),'',kind=pI64) /= 0_pI64) then gotCellData = .true. - startPos = endPos + 2_pI64 do while (index(fileContent(startPos:endPos),'',kind=pI64) == 0_pI64) - endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64 if(index(fileContent(startPos:endPos),'',kind=pI64) /= 0_pI64) then gotCoordinates = .true. @@ -415,10 +413,10 @@ subroutine readVTR(grid,geomSize,origin,microstructure) end do - if(.not. allocated(microstructure)) call IO_error(error_ID = 844, ext_msg='materialpoint not found') - if(size(microstructure) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(materialpoint)') - 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(.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found') + 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(grid<1)) call IO_error(error_ID = 844, ext_msg='grid') contains diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index f0d65a4a6..27dc5c1bd 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -22,42 +22,38 @@ module grid_damage_spectral implicit none private - type, private :: tNumerics + type :: tNumerics integer :: & - itmax !< max number of iterations + itmax !< maximum number of iterations real(pReal) :: & residualStiffness, & !< non-zero residual damage eps_damage_atol, & !< absolute tolerance for damage evolution eps_damage_rtol !< relative tolerance for damage evolution end type tNumerics - type(tNumerics), private :: num -!-------------------------------------------------------------------------------------------------- -! derived types - type(tSolutionParams), private :: params + type(tNumerics) :: num + type(tSolutionParams) :: params !-------------------------------------------------------------------------------------------------- ! PETSc data - SNES, private :: damage_snes - Vec, private :: solution_vec - PetscInt, private :: xstart, xend, ystart, yend, zstart, zend - real(pReal), private, dimension(:,:,:), allocatable :: & + SNES :: damage_snes + Vec :: solution_vec + PetscInt :: xstart, xend, ystart, yend, zstart, zend + real(pReal), dimension(:,:,:), allocatable :: & phi_current, & !< field of current damage phi_lastInc, & !< field of previous damage phi_stagInc !< field of staggered damage !-------------------------------------------------------------------------------------------------- ! reference diffusion tensor, mobility etc. - integer, private :: totalIter = 0 !< total iteration in current increment - real(pReal), dimension(3,3), private :: K_ref - real(pReal), private :: mu_ref + integer :: totalIter = 0 !< total iteration in current increment + real(pReal), dimension(3,3) :: K_ref + real(pReal) :: mu_ref public :: & grid_damage_spectral_init, & grid_damage_spectral_solution, & grid_damage_spectral_forward - private :: & - formResidual contains @@ -77,10 +73,10 @@ subroutine grid_damage_spectral_init character(len=pStringLen) :: & 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' - write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' + print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' + print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' !------------------------------------------------------------------------------------------------- ! 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) call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr) -!-------------------------------------------------------------------------------------------------- -! damage reference diffusion update call updateReference 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 VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr) if (solution%converged) & - write(6,'(/,a)') ' ... nonlocal damage converged .....................................' - write(6,'(/,a,f8.6,2x,f8.6,2x,e11.4,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& + print'(/,a)', ' ... nonlocal damage converged .....................................' + write(IO_STDOUT,'(/,a,f8.6,2x,f8.6,2x,e11.4,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& phi_min, phi_max, stagNorm - write(6,'(/,a)') ' ===========================================================================' - flush(6) + print'(/,a)', ' ===========================================================================' + flush(IO_STDOUT) end function grid_damage_spectral_solution diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 748692b3c..e1ceb42b0 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -122,7 +122,7 @@ subroutine grid_mech_FEM_init PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc - write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>'; flush(6) + print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT) !----------------------------------------------------------------------------------------------- ! debugging options @@ -130,13 +130,12 @@ subroutine grid_mech_FEM_init 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%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_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%itmin = num_grid%get_asInt ('itmin',defaultVal=1) num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) @@ -225,7 +224,7 @@ subroutine grid_mech_FEM_init !-------------------------------------------------------------------------------------------------- ! init fields 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' fileHandle = HDF5_openFile(fileName) @@ -254,7 +253,7 @@ subroutine grid_mech_FEM_init CHKERRQ(ierr) 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_volAvgLastInc,'C_volAvgLastInc') @@ -304,11 +303,11 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation !-------------------------------------------------------------------------------------------------- ! 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 - call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr) + call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) solution%converged = reason > 0 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_aim_lastInc = F_aim - !-------------------------------------------------------------------------------------------------- + !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F 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_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' fileHandle = HDF5_openFile(fileName,'w') @@ -476,13 +475,13 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i !-------------------------------------------------------------------------------------------------- ! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + print'(1/,a)', ' ... reporting .............................................................' + print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & 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,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + print'(/,a)', ' ===========================================================================' + flush(IO_STDOUT) end subroutine converged @@ -516,13 +515,13 @@ subroutine formResidual(da_local,x_local, & ! begin of new iteration newIteration: if (totalIter <= PETScIter) then 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) & - 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.)) - 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) - flush(6) + flush(IO_STDOUT) endif newIteration !-------------------------------------------------------------------------------------------------- @@ -541,7 +540,7 @@ subroutine formResidual(da_local,x_local, & !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(P_current,& + call utilities_constitutiveResponse(P_current,& P_av,C_volAvg,devNull, & F,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 9d301b30f..e7fb9bd19 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -42,8 +42,7 @@ module grid_mech_spectral_basic type(tNumerics) :: num ! numerics parameters. Better name? - logical, private:: & - debugRotation + logical, private :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -110,13 +109,13 @@ subroutine grid_mech_spectral_basic_init character(len=pStringLen) :: & 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' - write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012' + print*, 'Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' + 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' - write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' + print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' + print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' !------------------------------------------------------------------------------------------------- ! 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_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%itmin = num_grid%get_asInt ('itmin',defaultVal=1) 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 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' 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 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_volAvgLastInc,'C_volAvgLastInc') @@ -377,7 +375,7 @@ subroutine grid_mech_spectral_basic_restartWrite 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' fileHandle = HDF5_openFile(fileName,'w') @@ -437,13 +435,13 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm !-------------------------------------------------------------------------------------------------- ! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + print'(1/,a)', ' ... reporting .............................................................' + print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & 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,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + print'(/,a)', ' ===========================================================================' + flush(IO_STDOUT) end subroutine converged @@ -475,13 +473,13 @@ subroutine formResidual(in, F, & ! begin of new iteration newIteration: if (totalIter <= PETScIter) then 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) & - 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.)) - 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) - flush(6) + flush(IO_STDOUT) endif newIteration !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 5414d02a8..5cc7d1602 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -123,10 +123,10 @@ subroutine grid_mech_spectral_polarisation_init character(len=pStringLen) :: & 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' - write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' + print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' + print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' !------------------------------------------------------------------------------------------------ ! debugging options @@ -134,9 +134,8 @@ subroutine grid_mech_spectral_polarisation_init debugRotation = debug_grid%contains('rotation') !------------------------------------------------------------------------------------------------- -! read numerical parameters +! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - 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_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,:,:,:) 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' 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 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_volAvgLastInc,'C_volAvgLastInc') @@ -434,7 +433,7 @@ subroutine grid_mech_spectral_polarisation_restartWrite F => FandF_tau(0: 8,:,:,:) 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' fileHandle = HDF5_openFile(fileName,'w') @@ -498,15 +497,15 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm !-------------------------------------------------------------------------------------------------- ! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + print'(1/,a)', ' ... reporting .............................................................' + print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & 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,')' - 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,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + print'(/,a)', ' ===========================================================================' + flush(IO_STDOUT) end subroutine converged @@ -558,13 +557,13 @@ subroutine formResidual(in, FandF_tau, & ! begin of new iteration newIteration: if (totalIter <= PETScIter) then 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) & - 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.)) - 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) - flush(6) + flush(IO_STDOUT) endif newIteration !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 5456f40eb..49be5ad7e 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -23,20 +23,17 @@ module grid_thermal_spectral implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types - type(tSolutionParams) :: params - type :: tNumerics integer :: & - itmax !< maximum number of iterations + itmax !< maximum number of iterations real(pReal) :: & - eps_thermal_atol, & !< absolute tolerance for thermal equilibrium - eps_thermal_rtol !< relative tolerance for thermal equilibrium + eps_thermal_atol, & !< absolute tolerance for thermal equilibrium + eps_thermal_rtol !< relative tolerance for thermal equilibrium end type tNumerics type(tNumerics) :: num + type(tSolutionParams) :: params !-------------------------------------------------------------------------------------------------- ! PETSc data SNES :: thermal_snes @@ -74,13 +71,13 @@ subroutine grid_thermal_spectral_init class(tNode), pointer :: & 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' - write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' + print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' + 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%itmax = num_grid%get_asInt ('itmax', defaultVal=250) 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 call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,& - num_grid%get_asString('petsc_options',defaultVal=''),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -110,7 +106,7 @@ subroutine grid_thermal_spectral_init DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point grid(1),grid(2),grid(3), & ! global grid 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 thermal_grid,ierr) ! handle, error CHKERRQ(ierr) @@ -159,8 +155,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) timeinc_old !< increment in time of last increment integer :: i, j, k, cell type(tSolutionState) :: solution - class(tNode), pointer :: & - num_grid PetscInt :: devNull 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 VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr) if (solution%converged) & - write(6,'(/,a)') ' ... thermal conduction converged ..................................' - write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',& + print'(/,a)', ' ... thermal conduction converged ..................................' + write(IO_STDOUT,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',& T_min, T_max, stagNorm - write(6,'(/,a)') ' ===========================================================================' - flush(6) + print'(/,a)', ' ===========================================================================' + flush(IO_STDOUT) end function grid_thermal_spectral_solution diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 1dc66a23f..dc6430c72 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -208,10 +208,10 @@ subroutine spectral_utilities_init debugPETSc = debug_grid%contains('petsc') - if(debugPETSc) write(6,'(3(/,a),/)') & + if(debugPETSc) print'(3(/,a),/)', & ' Initializing PETSc with debug options: ', & 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) @@ -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' 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 @@ -506,8 +506,8 @@ subroutine utilities_fourierGammaConvolution(fieldAim) logical :: err - write(6,'(/,a)') ' ... doing gamma convolution ...............................................' - flush(6) + print'(/,a)', ' ... doing gamma convolution ...............................................' + flush(IO_STDOUT) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation (mechanical equilibrium) @@ -576,8 +576,8 @@ real(pReal) function utilities_divergenceRMS() integer :: i, j, k, ierr complex(pReal), dimension(3) :: rescaledGeom - write(6,'(/,a)') ' ... calculating divergence ................................................' - flush(6) + print'(/,a)', ' ... calculating divergence ................................................' + flush(IO_STDOUT) 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) :: rescaledGeom - write(6,'(/,a)') ' ... calculating curl ......................................................' - flush(6) + print'(/,a)', ' ... calculating curl ......................................................' + flush(IO_STDOUT) 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)) if(debugGeneral) then - write(6,'(/,a)') ' ... updating masked compliance ............................................' - write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& + print'(/,a)', ' ... updating masked compliance ............................................' + write(IO_STDOUT,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& transpose(temp99_Real)*1.0e-9_pReal - flush(6) + flush(IO_STDOUT) endif 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 write(formatString, '(i2)') size_reduced 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)) - 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') endif 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) 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 - flush(6) + flush(IO_STDOUT) endif 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), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF - write(6,'(/,a)') ' ... evaluating constitutive response ......................................' - flush(6) + print'(/,a)', ' ... evaluating constitutive response ......................................' + flush(IO_STDOUT) 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 call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) 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 if(present(rotation_BC)) & 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 - flush(6) + flush(IO_STDOUT) dPdF_max = 0.0_pReal dPdF_norm_max = 0.0_pReal @@ -1095,7 +1095,7 @@ subroutine utilities_saveReferenceStiffness fileUnit,ierr 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',& status='replace',access='stream',action='write',iostat=ierr) if(ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref') diff --git a/src/homogenization.f90 b/src/homogenization.f90 index cc8df77f4..0e7f1bf3a 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -186,7 +186,7 @@ subroutine homogenization_init materialpoint_F = materialpoint_F0 ! initialize to identity 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%subStepMinHomog = num_homogGeneric%get_asFloat('subStepMin', defaultVal=1.0e-3_pReal) diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index cdd7c05dd..f0485b244 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -95,7 +95,7 @@ module subroutine mech_RGC_init(num_homogMech) print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>' 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*, '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) enddo print*,' ' - flush(6) + flush(IO_STDOUT) endif #endif enddo @@ -376,7 +376,7 @@ module procedure mech_RGC_updateState '@ grain ',stresLoc(3),' in component ',stresLoc(1),stresLoc(2) print'(a,e15.8,a,i3,a,i2)',' Max residual: ',residMax, & ' @ iface ',residLoc(1),' in direction ',residLoc(2) - flush(6) + flush(IO_STDOUT) endif #endif @@ -388,7 +388,7 @@ module procedure mech_RGC_updateState mech_RGC_updateState = .true. #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) & - print*, '... done and happy'; flush(6) + print*, '... done and happy'; flush(IO_STDOUT) #endif !-------------------------------------------------------------------------------------------------- @@ -416,7 +416,7 @@ module procedure mech_RGC_updateState print'(a,e15.8,/)', ' Volume discrepancy: ', dst%volumeDiscrepancy(of) print'(a,e15.8)', ' Maximum relaxation rate: ', dst%relaxationRate_max(of) print'(a,e15.8,/)', ' Average relaxation rate: ', dst%relaxationRate_avg(of) - flush(6) + flush(IO_STDOUT) endif #endif @@ -429,7 +429,7 @@ module procedure mech_RGC_updateState #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) & - print'(a,/)', ' ... broken'; flush(6) + print'(a,/)', ' ... broken'; flush(IO_STDOUT) #endif return @@ -437,7 +437,7 @@ module procedure mech_RGC_updateState else ! proceed with computing the Jacobian and state update #ifdef DEBUG 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 @@ -499,7 +499,7 @@ module procedure mech_RGC_updateState print'(1x,100(e11.4,1x))',(smatrix(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' - flush(6) + flush(IO_STDOUT) endif #endif @@ -559,7 +559,7 @@ module procedure mech_RGC_updateState print'(1x,100(e11.4,1x))',(pmatrix(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' - flush(6) + flush(IO_STDOUT) endif #endif @@ -578,7 +578,7 @@ module procedure mech_RGC_updateState print'(1x,100(e11.4,1x))',(rmatrix(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' - flush(6) + flush(IO_STDOUT) endif #endif @@ -593,7 +593,7 @@ module procedure mech_RGC_updateState print'(1x,100(e11.4,1x))',(jmatrix(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' - flush(6) + flush(IO_STDOUT) endif #endif @@ -609,7 +609,7 @@ module procedure mech_RGC_updateState print'(1x,100(e11.4,1x))',(jnverse(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' - flush(6) + flush(IO_STDOUT) endif #endif @@ -625,7 +625,7 @@ module procedure mech_RGC_updateState !$OMP CRITICAL (write2out) 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)) - flush(6) + flush(IO_STDOUT) !$OMP END CRITICAL (write2out) endif @@ -636,7 +636,7 @@ module procedure mech_RGC_updateState print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of) enddo print*,' ' - flush(6) + flush(IO_STDOUT) endif #endif diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index 6c5f50b99..5138afa73 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -40,7 +40,7 @@ module subroutine mech_isostrain_init print'(/,a)', ' <<<+- homogenization_mech_isostrain init -+>>>' 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 diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90 index d9426ef50..3cbec5911 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mech_none.f90 @@ -21,7 +21,7 @@ module subroutine mech_none_init print'(/,a)', ' <<<+- homogenization_mech_none init -+>>>' 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) if (homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index d8f25f8b8..23f348831 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -49,7 +49,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin myKinematics = kinematics_active('cleavage_opening',kinematics_length) Ninstance = count(myKinematics) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return phases => config_material%get('phase') diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 3b04e37c1..660483b90 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -52,7 +52,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi myKinematics = kinematics_active('slipplane_opening',kinematics_length) Ninstance = count(myKinematics) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return phases => config_material%get('phase') diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 652713aa4..93a48e035 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -42,7 +42,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi myKinematics = kinematics_active('thermal_expansion',kinematics_length) Ninstance = count(myKinematics) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return phases => config_material%get('phase') diff --git a/src/lattice.f90 b/src/lattice.f90 index d0ac07ed6..69305c839 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -457,7 +457,7 @@ subroutine lattice_init phase, & elasticity - print'(/,a)', ' <<<+- lattice init -+>>>'; flush(6) + print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) phases => config_material%get('phase') Nphases = phases%length diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index a117a6d88..9696572e5 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -52,7 +52,7 @@ subroutine discretization_marc_init type(tElement) :: elem integer, dimension(:), allocatable :: & - microstructureAt + materialAt integer:: & Nnodes, & !< total number of nodes in the mesh Nelems, & !< total number of elements in the mesh @@ -70,7 +70,7 @@ subroutine discretization_marc_init class(tNode), pointer :: & num_commercialFEM - write(6,'(/,a)') ' <<<+- discretization_marc init -+>>>'; flush(6) + print'(/,a)', ' <<<+- discretization_marc init -+>>>'; flush(6) !--------------------------------------------------------------------------------- ! 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 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) 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,& elem%nIPs*nElems]),node0_cell) - call discretization_init(microstructureAt,& + call discretization_init(materialAt,& IP_reshaped,& node0_cell) @@ -172,7 +172,7 @@ end subroutine writeGeometry !-------------------------------------------------------------------------------------------------- !> @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 real(pReal), dimension(:,:), allocatable, intent(out) :: & @@ -180,7 +180,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt) integer, dimension(:,:), allocatable, intent(out) :: & connectivity_elem integer, dimension(:), allocatable, intent(out) :: & - microstructureAt + materialAt integer :: & fileFormatVersion, & @@ -226,9 +226,9 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt) connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,inputFile) - call inputRead_microstructure(microstructureAt, & - nElems,elem%nNodes,nameElemSet,mapElemSet,& - initialcondTableStyle,inputFile) + call inputRead_material(materialAt, & + nElems,elem%nNodes,nameElemSet,mapElemSet,& + initialcondTableStyle,inputFile) end subroutine inputRead @@ -675,13 +675,13 @@ end function inputRead_connectivityElem !-------------------------------------------------------------------------------------------------- -!> @brief Store microstructure ID +!> @brief Store material ID !-------------------------------------------------------------------------------------------------- -subroutine inputRead_microstructure(microstructureAt,& - nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileContent) +subroutine inputRead_material(materialAt,& + nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileContent) integer, dimension(:), allocatable, intent(out) :: & - microstructureAt + materialAt integer, intent(in) :: & nElem, & 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 - allocate(microstructureAt(nElem),source=0) + allocate(materialAt(nElem),source=0) do l = 1, size(fileContent) 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 do i = 1,contInts(1) e = mesh_FEM2DAMASK_elem(contInts(1+i)) - microstructureAt(e) = myVal + materialAt(e) = myVal enddo if (initialcondTableStyle == 0) m = m + 1 enddo @@ -723,9 +723,9 @@ subroutine inputRead_microstructure(microstructureAt,& endif 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), & nodePos(1:3,3) - nodePos(1:3,1)) case (4) ! 3D 8node - ! for this cell type we get the normal of the quadrilateral face as an average of - ! four normals of triangular subfaces; since the face consists only of two triangles, - ! the sum has to be divided by two; this whole prcedure tries to compensate for - ! probable non-planar cell surfaces + ! Get the normal of the quadrilateral face as the average of four normals of triangular + ! subfaces. Since the face consists only of two triangles, the sum has to be divided + ! by two. This procedure tries to compensate for probable non-planar cell surfaces IPareaNormal(1:3,f,i,e) = 0.0_pReal do n = 1, m IPareaNormal(1:3,f,i,e) = IPareaNormal(1:3,f,i,e) & diff --git a/src/material.f90 b/src/material.f90 index 835de2fc1..6688a0ae5 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -164,7 +164,7 @@ subroutine material_init(restart) material_homogenization character(len=pStringLen) :: sectionName - print'(/,a)', ' <<<+- material init -+>>>'; flush(6) + print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) phases => config_material%get('phase') allocate(material_name_phase(phases%length)) @@ -181,10 +181,10 @@ subroutine material_init(restart) enddo call material_parseMicrostructure - print*, ' Microstructure parsed' + print*, 'Microstructure parsed' call material_parseHomogenization - print*, ' Homogenization parsed' + print*, 'Homogenization parsed' 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 +! ToDo: This should be done in homogenization !-------------------------------------------------------------------------------------------------- subroutine material_parseHomogenization @@ -320,100 +321,78 @@ end subroutine material_parseHomogenization !-------------------------------------------------------------------------------------------------- subroutine material_parseMicrostructure - class(tNode), pointer :: microstructure, & !> pointer to microstructure list - constituentsInMicrostructure, & !> pointer to a microstructure list item - constituents, & !> pointer to constituents list - constituent, & !> pointer to each constituent + class(tNode), pointer :: microstructures, & !> list of microstructures + microstructure, & !> microstructure definition + constituents, & !> list of constituents + constituent, & !> constituent definition phases, & homogenization integer, dimension(:), allocatable :: & - CounterPhase, & - CounterHomogenization - - - real(pReal), dimension(:,:), allocatable :: & - microstructure_fraction !< vol fraction of each constituent in microstrcuture + counterPhase, & + counterHomogenization + real(pReal) :: & + frac integer :: & e, & i, & m, & 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') - phases => config_material%get('phase') - microstructure => config_material%get('microstructure') - allocate(microstructure_Nconstituents(microstructure%length), source = 0) - - 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') + allocate(microstructure_Nconstituents(microstructures%length),source=0) + do m = 1, microstructures%length + microstructure => microstructures%get(m) + constituents => microstructure%get('constituents') microstructure_Nconstituents(m) = constituents%length enddo - - microstructure_maxNconstituents = maxval(microstructure_Nconstituents) - allocate(microstructure_fraction(microstructure_maxNconstituents,microstructure%length), source =0.0_pReal) - allocate(material_phaseAt(microstructure_maxNconstituents,discretization_nElem), source =0) - allocate(material_orientation0(microstructure_maxNconstituents,discretization_nIP,discretization_nElem)) - allocate(material_homogenizationAt(discretization_nElem)) + maxNconstituents = maxval(microstructure_Nconstituents) + + allocate(material_homogenizationAt(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(CounterHomogenization(homogenization%length),source=0) + allocate(material_orientation0(maxNconstituents,discretization_nIP,discretization_nElem)) + + 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 - constituentsInMicrostructure => microstructure%get(m) - constituents => constituentsInMicrostructure%get('constituents') + do e = 1, discretization_nElem + microstructure => microstructures%get(discretization_microstructureAt(e)) + 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 constituent => constituents%get(c) - microstructure_fraction(c,m) = constituent%get_asFloat('fraction') - enddo - if (dNeq(sum(microstructure_fraction(:,m)),1.0_pReal)) call IO_error(153,ext_msg='constituent') - 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 - 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) + frac = frac + constituent%get_asFloat('fraction') + + material_phaseAt(c,e) = phases%getIndex(constituent%get_asString('phase')) + do i = 1, discretization_nIP + counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1 + material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e)) + + call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('orientation',requiredSize=4)) enddo + enddo + if (dNeq(frac,1.0_pReal)) call IO_error(153,ext_msg='constituent') + 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 diff --git a/src/math.f90 b/src/math.f90 index b835a35b2..163f4df6a 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -91,7 +91,7 @@ subroutine math_init class(tNode), pointer :: & num_generic - print'(/,a)', ' <<<+- math init -+>>>'; flush(6) + print'(/,a)', ' <<<+- math init -+>>>'; flush(IO_STDOUT) num_generic => config_numerics%get('generic',defaultVal=emptyDict) randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 2fefe73a9..bfa8d22ce 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -78,7 +78,7 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) 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 @@ -208,30 +208,30 @@ program DAMASK_mesh errorID = 0 checkLoadcases: do currentLoadCase = 1, size(loadCases) write (loadcase_string, '(i0)' ) currentLoadCase - write(6,'(1x,a,i6)') 'load case: ', currentLoadCase + print'(a,i0)', ' load case: ', currentLoadCase if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & - write(6,'(2x,a)') 'drop guessing along trajectory' + print'(a)', ' drop guessing along trajectory' do field = 1, nActiveFields select case (loadCases(currentLoadCase)%fieldBC(field)%ID) case(FIELD_MECH_ID) - write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label) + print'(a)', ' Field '//trim(FIELD_MECH_label) end select do faceSet = 1, mesh_Nboundaries do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents 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, & ' Value ', loadCases(currentLoadCase)%fieldBC(field)% & componentBC(component)%Value(faceSet) 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 - 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 - write(6,'(2x,a,i5)') 'output frequency: ', & + print'(a,i5)', ' output frequency: ', & loadCases(currentLoadCase)%outputfrequency if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message enddo checkLoadcases @@ -290,8 +290,8 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! report begin of new step - write(6,'(/,a)') ' ###########################################################################' - write(6,'(1x,a,es12.5,6(a,i0))')& + print'(/,a)', ' ###########################################################################' + print'(1x,a,es12.5,6(a,i0))',& 'Time', time, & 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& '-', stepFraction, '/', subStepFactor**cutBackLevel,& @@ -299,7 +299,7 @@ program DAMASK_mesh write(incInfo,'(4(a,i0))') & 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& '-',stepFraction, '/', subStepFactor**cutBackLevel - flush(6) + flush(IO_STDOUT) !-------------------------------------------------------------------------------------------------- ! forward fields @@ -338,7 +338,7 @@ program DAMASK_mesh cutBack = .False. if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if (cutBackLevel < maxCutBack) then ! do cut back - write(6,'(/,a)') ' cut back detected' + print'(/,a)', ' cut back detected' cutBack = .True. stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator cutBackLevel = cutBackLevel + 1 @@ -360,13 +360,13 @@ program DAMASK_mesh cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc if (all(solres(:)%converged)) then - write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' converged' + print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged' else - write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' NOT converged' - endif; flush(6) + print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged' + endif; flush(IO_STDOUT) 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) endif @@ -378,7 +378,7 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation - write(6,'(/,a)') ' ###########################################################################' + print'(/,a)', ' ###########################################################################' if (worldrank == 0) close(statUnit) call quit(0) ! no complains ;) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 65b0bc6fe..4cc5e5468 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -37,7 +37,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine FEM_quadrature_init - write(6,'(/,a)') ' <<<+- FEM_quadrature init -+>>>'; flush(6) + print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6) !-------------------------------------------------------------------------------------------------- ! 2D linear diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 551a863b6..4d9786112 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -110,7 +110,7 @@ subroutine FEM_utilities_init PetscErrorCode :: ierr - write(6,'(/,a)') ' <<<+- FEM_utilities init -+>>>' + print'(/,a)', ' <<<+- FEM_utilities init -+>>>' num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2) @@ -118,11 +118,11 @@ subroutine FEM_utilities_init debug_mesh => config_debug%get('mesh',defaultVal=emptyList) debugPETSc = debug_mesh%contains('petsc') - if(debugPETSc) write(6,'(3(/,a),/)') & + if(debugPETSc) print'(3(/,a),/)', & ' Initializing PETSc with debug options: ', & trim(PETScDebug), & ' add more using the PETSc_Options keyword in numerics.yaml ' - flush(6) + flush(IO_STDOUT) call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr) CHKERRQ(ierr) if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) @@ -158,7 +158,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) PetscErrorCode :: ierr - write(6,'(/,a)') ' ... evaluating constitutive response ......................................' + print'(/,a)', ' ... evaluating constitutive response ......................................' call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 54854ce68..b57a4b332 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -83,7 +83,7 @@ subroutine discretization_mesh_init(restart) num_mesh integer :: integrationOrder !< order of quadrature rule required - write(6,'(/,a)') ' <<<+- discretization_mesh init -+>>>' + print'(/,a)', ' <<<+- discretization_mesh init -+>>>' !-------------------------------------------------------------------------------- ! read numerics parameter diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 546897c5b..de1f0c687 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -110,7 +110,7 @@ subroutine FEM_mech_init(fieldBC) class(tNode), pointer :: & 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 @@ -318,8 +318,8 @@ type(tSolutionState) function FEM_mech_solution( & CHKERRQ(ierr) endif - write(6,'(/,a)') ' ===========================================================================' - flush(6) + print'(/,a)', ' ===========================================================================' + flush(IO_STDOUT) 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) CHKERRQ(ierr) 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 = ', & 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 - flush(6) + flush(IO_STDOUT) end subroutine FEM_mech_converged diff --git a/src/parallelization.f90 b/src/parallelization.f90 index 341b620c9..3bb1f6af5 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -3,8 +3,8 @@ !> @brief Inquires variables related to parallelization (openMP, MPI) !-------------------------------------------------------------------------------------------------- module parallelization - use prec - use, intrinsic :: iso_fortran_env + use, intrinsic :: ISO_fortran_env, only: & + OUTPUT_UNIT #ifdef PETSc #include @@ -12,6 +12,8 @@ module parallelization #endif !$ use OMP_LIB + use prec + implicit none private @@ -29,14 +31,14 @@ contains !-------------------------------------------------------------------------------------------------- subroutine parallelization_init - integer :: err, typeSize + integer :: err, typeSize !$ integer :: got_env, DAMASK_NUM_THREADS, threadLevel !$ character(len=6) NumThreadsString #ifdef PETSc PetscErrorCode :: petsc_err #else - print'(/,a)', ' <<<+- parallelization init -+>>>'; flush(6) + print'(/,a)', ' <<<+- parallelization init -+>>>'; flush(OUTPUT_UNIT) #endif #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' #endif - mainProcess: if (worldrank == 0) then - if (output_unit /= 6) error stop 'STDOUT != 6' - if (error_unit /= 0) error stop 'STDERR != 0' - else mainProcess - 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 - + if (worldrank /= 0) then + close(OUTPUT_UNIT) ! disable output + open(OUTPUT_UNIT,file='/dev/null',status='replace') ! close() alone will leave some temp files in cwd + endif !$ call get_environment_variable(name='DAMASK_NUM_THREADS',value=NumThreadsString,STATUS=got_env) !$ if(got_env /= 0) then diff --git a/src/prec.f90 b/src/prec.f90 index c8700089b..cec4928ba 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -8,7 +8,7 @@ !-------------------------------------------------------------------------------------------------- module prec use, intrinsic :: IEEE_arithmetic - use, intrinsic :: ISO_C_Binding + use, intrinsic :: ISO_C_binding implicit none public diff --git a/src/results.f90 b/src/results.f90 index 686183919..aec90d7be 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -65,7 +65,7 @@ subroutine results_init(restart) 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*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL diff --git a/src/rotations.f90 b/src/rotations.f90 index 72046a965..e19513866 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -104,7 +104,7 @@ contains subroutine rotations_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*, 'https://doi.org/10.1088/0965-0393/23/8/083501' diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 6dd58fe5b..7911d6d0a 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -53,7 +53,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) mySources = source_active('damage_anisoBrittle',source_length) Ninstance = count(mySources) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return phases => config_material%get('phase') diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 37681a23f..52189c839 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -47,7 +47,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources) mySources = source_active('damage_anisoDuctile',source_length) Ninstance = count(mySources) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return phases => config_material%get('phase') diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index c6d0ada99..714e71ef1 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -43,7 +43,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources) mySources = source_active('damage_isoBrittle',source_length) Ninstance = count(mySources) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return phases => config_material%get('phase') diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 1c1c53fd0..493183d75 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -45,7 +45,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) mySources = source_active('damage_isoDuctile',source_length) Ninstance = count(mySources) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return phases => config_material%get('phase') diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 60fde7f28..5cc740424 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -41,7 +41,7 @@ module function source_thermal_dissipation_init(source_length) result(mySources) mySources = source_active('thermal_dissipation',source_length) Ninstance = count(mySources) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return phases => config_material%get('phase') diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index e7bfea254..4a644f53b 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -45,7 +45,7 @@ module function source_thermal_externalheat_init(source_length) result(mySources mySources = source_active('thermal_externalheat',source_length) Ninstance = count(mySources) - print'(a,i2)', ' # instances: ',Ninstance; flush(6) + print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) if(Ninstance == 0) return phases => config_material%get('phase')