From 83ccd7c982d63b9c0ec45f43d4161dc825033d31 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Mar 2016 21:37:22 +0100 Subject: [PATCH 01/24] introduced (in)equality comparison for double precision --- code/prec.f90 | 30 ++++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/code/prec.f90 b/code/prec.f90 index 0f22cf649..129681fd5 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -172,9 +172,9 @@ end subroutine prec_init !-------------------------------------------------------------------------------------------------- !> @brief figures out if a floating point number is NaN -! basically just a small wrapper, because gfortran < 4.9 does not have the IEEE module +! basically just a small wrapper, because gfortran < 5.0 does not have the IEEE module !-------------------------------------------------------------------------------------------------- -logical elemental function prec_isNaN(a) +logical elemental pure function prec_isNaN(a) implicit none real(pReal), intent(in) :: a @@ -187,4 +187,30 @@ logical elemental function prec_isNaN(a) #endif end function prec_isNaN + + +!-------------------------------------------------------------------------------------------------- +!> @brief equality comparison for double precision +! replaces "==" but for certain (absolute) tolerance. Counterpart to dNeq +!-------------------------------------------------------------------------------------------------- +logical elemental pure function dEq(a,b,tol) + real(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal), parameter :: eps = 1.0e-15_pReal + dEq = merge(.True., .False.,abs(a-b) < merge(tol,eps,present(tol))) +end function dEq + + + +!-------------------------------------------------------------------------------------------------- +!> @brief inequality comparison for double precision +! replaces "!=" but for certain (absolute) tolerance. Counterpart to dEq +!-------------------------------------------------------------------------------------------------- +logical elemental pure function dNeq(a,b,tol) + real(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal), parameter :: eps = 1.0e-15_pReal + dNeq = merge(.True., .False.,abs(a-b) < merge(tol,eps,present(tol))) +end function dNeq + end module prec From 69c1cb71de23f0fb1fc5b616fddedc05c3e11c05 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 9 Mar 2016 04:12:09 +0100 Subject: [PATCH 02/24] updated version information after successful test of v2.0.0-6-gb5302ec --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 2bc4b2cfc..d9e8f873a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-3-g618e352 +v2.0.0-6-gb5302ec From bb9ee3910fac62010ab646cf3aaa841d0ccf6d7e Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 9 Mar 2016 16:12:54 +0100 Subject: [PATCH 03/24] updated version information after successful test of v2.0.0-7-g2296929 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 2bc4b2cfc..4c6293bf0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-3-g618e352 +v2.0.0-7-g2296929 From 167b2783276b23fdbd8526cbf4d74fc9012cd42b Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 12 Mar 2016 05:21:45 +0100 Subject: [PATCH 04/24] updated version information after successful test of v2.0.0-16-g344c6f6 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 99c96eda6..c02195acf 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-15-g3898fb2 +v2.0.0-16-g344c6f6 From 2dbdf9a72f28f488746b598b9914ae008b03b455 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Mar 2016 18:21:25 +0100 Subject: [PATCH 05/24] all symbols need to be representable in unicode --- lib/damask/util.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/damask/util.py b/lib/damask/util.py index 77bc30403..9ce42759a 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -104,7 +104,8 @@ class extendableOption(Option): class backgroundMessage(threading.Thread): """reporting with animation to indicate progress""" - choices = {'bounce': ['_','o','O','°','¯','¯','°','O','o',], + choices = {'bounce': ['_', 'o', 'O', u'\u00B0', + u'\u203e',u'\u203e',u'\u00B0','O','o','_'], 'circle': [u'\u25f4',u'\u25f5',u'\u25f6',u'\u25f7'], 'hexagon': [u'\u2b22',u'\u2b23'], 'square': [u'\u2596',u'\u2598',u'\u259d',u'\u2597'], From 80d2bdbfde1c55fbf3332eeedfd8565efa19de65 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 14 Mar 2016 04:12:59 +0100 Subject: [PATCH 06/24] updated version information after successful test of v2.0.0-18-g2dbdf9a --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index c02195acf..6a337695e 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-16-g344c6f6 +v2.0.0-18-g2dbdf9a From b40e7a26124a59b8d38240d4b49a2012e0a42f3d Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 14 Mar 2016 16:11:55 +0100 Subject: [PATCH 07/24] updated version information after successful test of v2.0.0-19-g80d2bdb --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 6a337695e..9f18a71da 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-18-g2dbdf9a +v2.0.0-19-g80d2bdb From 6197e839c942b5dc91fb7243bc0423fff39b3264 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 15 Mar 2016 14:52:10 -0400 Subject: [PATCH 08/24] fixed missing file extension when adding in-place --- processing/post/vtk_addRectilinearGridData.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/vtk_addRectilinearGridData.py b/processing/post/vtk_addRectilinearGridData.py index 922b4b696..3ca54d84e 100755 --- a/processing/post/vtk_addRectilinearGridData.py +++ b/processing/post/vtk_addRectilinearGridData.py @@ -151,7 +151,7 @@ for name in filenames: writer = vtk.vtkXMLRectilinearGridWriter() writer.SetDataModeToBinary() writer.SetCompressorTypeToZLib() - writer.SetFileName(os.path.splitext(options.vtk)[0]+('' if options.inplace else '_added.vtr')) + writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtr' if options.inplace else '_added.vtr')) if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid) else: writer.SetInputData(rGrid) writer.Write() From 9704eb6926f6be905128ae0552dec146fa811eb7 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 15 Mar 2016 17:23:01 -0400 Subject: [PATCH 09/24] small correction to variable declaration layout --- code/lattice.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/lattice.f90 b/code/lattice.f90 index f5717e8e1..05a123125 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -2182,7 +2182,7 @@ pure function lattice_qDisorientation(Q1, Q2, struct) real(pReal), dimension(4) :: lattice_qDisorientation real(pReal), dimension(4), intent(in) :: & Q1, & ! 1st orientation - Q2 ! 2nd orientation + Q2 ! 2nd orientation integer(kind(LATTICE_undefined_ID)), optional, intent(in) :: & ! if given, symmetries between the two orientation will be considered struct From 93466273ffd03db9cc03e6fadf02727ba3967d19 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 16 Mar 2016 04:13:16 +0100 Subject: [PATCH 10/24] updated version information after successful test of v2.0.0-31-g6197e83 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 9f18a71da..27bd816bb 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-19-g80d2bdb +v2.0.0-31-g6197e83 From ce26ad06da8cbc945444412d42220b4755562730 Mon Sep 17 00:00:00 2001 From: Aritra Chakraborty Date: Wed, 16 Mar 2016 14:55:56 -0400 Subject: [PATCH 11/24] polishing and introduction of locally derived grid in FFT subroutine --- processing/post/addCurl.py | 39 +++++++++++++------------------- processing/post/addDivergence.py | 36 ++++++++++++----------------- 2 files changed, 31 insertions(+), 44 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index ac9715c36..6789b27d5 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -10,40 +10,35 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def curlFFT(geomdim,field): + grid = np.array(np.shape(field)[2::-1]) N = grid.prod() # field size n = np.array(np.shape(field)[3:]).prod() # data size - if n == 3: - dataType = 'vector' - elif n == 9: - dataType = 'tensor' + if n == 3: dataType = 'vector' + elif n == 9: dataType = 'tensor' field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2)) curl_fourier = np.zeros(field_fourier.shape,'c16') # differentiation in Fourier space k_s = np.zeros([3],'i') - TWOPIIMG = (0.0+2.0j*math.pi) + TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if(grid[2]%2==0 and i == grid[2]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[0]=0 - elif (i > grid[2]//2): - k_s[0] = k_s[0] - grid[2] + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if(grid[1]%2==0 and j == grid[1]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[1]=0 - elif (j > grid[1]//2): - k_s[1] = k_s[1] - grid[1] + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if(grid[0]%2==0 and k == grid[0]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[2]=0 + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order - xi = np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16') if dataType == 'tensor': for l in xrange(3): curl_fourier[i,j,k,0,l] = ( field_fourier[i,j,k,l,2]*xi[1]\ @@ -100,10 +95,8 @@ if options.vector is None and options.tensor is None: if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name,buffered = False) - except: - continue + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------ @@ -161,9 +154,9 @@ for name in filenames: stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(curlFFT(size[::-1], # we need to reverse order here, because x - table.data[:,data['column'][i]:data['column'][i]+data['dim']]. # is fastest,ie rightmost, but leftmost in - reshape([grid[2],grid[1],grid[0]]+data['shape']))) # our x,y,z notation + stack.append(curlFFT(size[::-1], # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + table.data[:,data['column'][i]:data['column'][i]+data['dim']]. + reshape([grid[2],grid[1],grid[0]]+data['shape']))) # ------------------------------------------ output result ----------------------------------------- diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 75708a603..99b56ead2 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -10,6 +10,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def divFFT(geomdim,field): + grid = np.array(np.shape(field)[2::-1]) N = grid.prod() # field size n = np.array(np.shape(field)[3:]).prod() # data size @@ -18,27 +19,22 @@ def divFFT(geomdim,field): # differentiation in Fourier space k_s=np.zeros([3],'i') - TWOPIIMG = (0.0+2.0j*math.pi) + TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if(grid[2]%2==0 and i == grid[2]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[0]=0 - elif (i > grid[2]//2): - k_s[0] = k_s[0] - grid[2] + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if(grid[1]%2==0 and j == grid[1]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[1]=0 - elif (j > grid[1]//2): - k_s[1] = k_s[1] - grid[1] + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if(grid[0]%2==0 and k == grid[0]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[2]=0 + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - xi=np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16') + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order if n == 9: # tensor, 3x3 -> 3 for l in xrange(3): div_fourier[i,j,k,l] = sum(field_fourier[i,j,k,l,0:3]*xi) *TWOPIIMG @@ -85,10 +81,8 @@ if options.vector is None and options.tensor is None: if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name,buffered = False) - except: - continue + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------ @@ -140,16 +134,16 @@ for name in filenames: maxcorner = np.array(map(max,coords)) grid = np.array(map(len,coords),'i') size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) - size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other spacings # ------------------------------------------ process value field ----------------------------------- stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(divFFT(size[::-1], # we need to reverse order here, because x - table.data[:,data['column'][i]:data['column'][i]+data['dim']]. # is fastest,ie rightmost, but leftmost in - reshape([grid[2],grid[1],grid[0]]+data['shape']))) # our x,y,z notation + stack.append(divFFT(size[::-1], # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + table.data[:,data['column'][i]:data['column'][i]+data['dim']]. + reshape([grid[2],grid[1],grid[0]]+data['shape']))) # ------------------------------------------ output result ----------------------------------------- @@ -158,4 +152,4 @@ for name in filenames: # ------------------------------------------ output finalization ----------------------------------- - table.close() # close input ASCII table (works for stdin) \ No newline at end of file + table.close() # close input ASCII table (works for stdin) From 3a267dc7ffc4f65d18bd30f9743c171e9d389d41 Mon Sep 17 00:00:00 2001 From: Aritra Chakraborty Date: Wed, 16 Mar 2016 15:02:38 -0400 Subject: [PATCH 12/24] FFT-based gradient calculation of scalar or vector field data (sibling of addCurl and addDivergence) --- processing/post/addGradient.py | 157 +++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100755 processing/post/addGradient.py diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py new file mode 100755 index 000000000..0d7caa829 --- /dev/null +++ b/processing/post/addGradient.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python +# -*- coding: UTF-8 no BOM -*- + +import os,sys,math +import numpy as np +from optparse import OptionParser +import damask + +scriptName = os.path.splitext(os.path.basename(__file__))[0] +scriptID = ' '.join([scriptName,damask.version]) + +def gradFFT(geomdim,field): + grid = np.array(np.shape(field)[2::-1]) + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size + if n == 3: dataType = 'vector' + elif n == 1: dataType = 'scalar' + + field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2)) + grad_fourier = np.zeros(field_fourier.shape+(3,),'c16') + +# differentiation in Fourier space + k_s = np.zeros([3],'i') + TWOPIIMG = 2.0j*math.pi + for i in xrange(grid[2]): + k_s[0] = i + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif i > grid[2]//2: k_s[0] -= grid[2] + + for j in xrange(grid[1]): + k_s[1] = j + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif j > grid[1]//2: k_s[1] -= grid[1] + + for k in xrange(grid[0]//2+1): + k_s[2] = k + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field order + + grad_fourier[i,j,k,0,:] = field_fourier[i,j,k,0]*xi *TWOPIIMG # vector field from scalar data + + if dataType == 'vector': + grad_fourier[i,j,k,1,:] = field_fourier[i,j,k,1]*xi *TWOPIIMG # tensor field from vector data + grad_fourier[i,j,k,2,:] = field_fourier[i,j,k,2]*xi *TWOPIIMG + + return np.fft.fftpack.irfftn(grad_fourier,axes=(0,1,2)).reshape([N,3*n]) + + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- + +parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ +Add column(s) containing gradient of requested column(s). +Operates on periodic ordered three-dimensional data sets. +Deals with both vector- and scalar fields. + +""", version = scriptID) + +parser.add_option('-c','--coordinates', + dest = 'coords', + type = 'string', metavar='string', + help = 'column heading for coordinates [%default]') +parser.add_option('-v','--vector', + dest = 'vector', + action = 'extend', metavar = '', + help = 'heading of columns containing vector field values') +parser.add_option('-s','--scalar', + dest = 'scalar', + action = 'extend', metavar = '', + help = 'heading of columns containing scalar field values') + +parser.set_defaults(coords = 'ipinitialcoord', + ) + +(options,filenames) = parser.parse_args() + +if options.vector is None and options.scalar is None: + parser.error('no data column specified.') + +# --- loop over input files ------------------------------------------------------------------------- + +if filenames == []: filenames = [None] + +for name in filenames: + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue + damask.util.report(scriptName,name) + +# ------------------------------------------ read header ------------------------------------------ + + table.head_read() + +# ------------------------------------------ sanity checks ---------------------------------------- + + items = { + 'scalar': {'dim': 1, 'shape': [1], 'labels':options.scalar, 'active':[], 'column': []}, + 'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []}, + } + errors = [] + remarks = [] + column = {} + + if table.label_dimension(options.coords) != 3: errors.append('coordinates {} are not a vector.'.format(options.coords)) + else: colCoord = table.label_index(options.coords) + + for type, data in items.iteritems(): + for what in (data['labels'] if data['labels'] is not None else []): + dim = table.label_dimension(what) + if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type)) + else: + items[type]['active'].append(what) + items[type]['column'].append(table.label_index(what)) + + if remarks != []: damask.util.croak(remarks) + if errors != []: + damask.util.croak(errors) + table.close(dismiss = True) + continue + +# ------------------------------------------ assemble header -------------------------------------- + + table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) + for type, data in items.iteritems(): + for label in data['active']: + table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in xrange(3 * data['dim'])]) # extend ASCII header with new labels # grad increases the field dimension by one + table.head_write() + +# --------------- figure out size and grid --------------------------------------------------------- + + table.data_readArray() + + coords = [np.unique(table.data[:,colCoord+i]) for i in xrange(3)] + mincorner = np.array(map(min,coords)) + maxcorner = np.array(map(max,coords)) + grid = np.array(map(len,coords),'i') + size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) + +# ------------------------------------------ process value field ----------------------------------- + + stack = [table.data] + for type, data in items.iteritems(): + for i,label in enumerate(data['active']): + stack.append(gradFFT(size[::-1], # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + table.data[:,data['column'][i]:data['column'][i]+data['dim']]. + reshape([grid[2],grid[1],grid[0]]+data['shape']))) + +# ------------------------------------------ output result ----------------------------------------- + + if len(stack) > 1: table.data = np.hstack(tuple(stack)) + table.data_writeArray('%.12g') + +# ------------------------------------------ output finalization ----------------------------------- + + table.close() # close input ASCII table (works for stdin) From 36c703ecd34f0c9033f372e6799dd831b6daede7 Mon Sep 17 00:00:00 2001 From: Aritra Chakraborty Date: Wed, 16 Mar 2016 15:12:53 -0400 Subject: [PATCH 13/24] shortened excessively long lines --- processing/post/addCurl.py | 15 ++++++++------- processing/post/addDivergence.py | 25 +++++++++++++------------ processing/post/addGradient.py | 23 ++++++++++++----------- 3 files changed, 33 insertions(+), 30 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 6789b27d5..d8b1ee025 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -11,8 +11,8 @@ scriptID = ' '.join([scriptName,damask.version]) def curlFFT(geomdim,field): grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - n = np.array(np.shape(field)[3:]).prod() # data size + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size if n == 3: dataType = 'vector' elif n == 9: dataType = 'tensor' @@ -25,19 +25,19 @@ def curlFFT(geomdim,field): TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order if dataType == 'tensor': for l in xrange(3): @@ -154,7 +154,8 @@ for name in filenames: stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(curlFFT(size[::-1], # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + stack.append(curlFFT(size[::-1], table.data[:,data['column'][i]:data['column'][i]+data['dim']]. reshape([grid[2],grid[1],grid[0]]+data['shape']))) diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 99b56ead2..aadaceabf 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -11,34 +11,34 @@ scriptID = ' '.join([scriptName,damask.version]) def divFFT(geomdim,field): grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - n = np.array(np.shape(field)[3:]).prod() # data size + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2)) - div_fourier = np.zeros(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector + div_fourier = np.zeros(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector # differentiation in Fourier space k_s=np.zeros([3],'i') TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order - if n == 9: # tensor, 3x3 -> 3 + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order + if n == 9: # tensor, 3x3 -> 3 for l in xrange(3): div_fourier[i,j,k,l] = sum(field_fourier[i,j,k,l,0:3]*xi) *TWOPIIMG - elif n == 3: # vector, 3 -> 1 + elif n == 3: # vector, 3 -> 1 div_fourier[i,j,k] = sum(field_fourier[i,j,k,0:3]*xi) *TWOPIIMG return np.fft.fftpack.irfftn(div_fourier,axes=(0,1,2)).reshape([N,n/3]) @@ -76,7 +76,7 @@ parser.set_defaults(coords = 'ipinitialcoord', if options.vector is None and options.tensor is None: parser.error('no data column specified.') -# --- loop over input files ------------------------------------------------------------------------- +# --- loop over input files ------------------------------------------------------------------------ if filenames == []: filenames = [None] @@ -134,14 +134,15 @@ for name in filenames: maxcorner = np.array(map(max,coords)) grid = np.array(map(len,coords),'i') size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) - size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other spacings + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones # ------------------------------------------ process value field ----------------------------------- stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(divFFT(size[::-1], # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + stack.append(divFFT(size[::-1], table.data[:,data['column'][i]:data['column'][i]+data['dim']]. reshape([grid[2],grid[1],grid[0]]+data['shape']))) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 0d7caa829..555e587de 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -11,8 +11,8 @@ scriptID = ' '.join([scriptName,damask.version]) def gradFFT(geomdim,field): grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - n = np.array(np.shape(field)[3:]).prod() # data size + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size if n == 3: dataType = 'vector' elif n == 1: dataType = 'scalar' @@ -24,24 +24,24 @@ def gradFFT(geomdim,field): TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field order + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field order - grad_fourier[i,j,k,0,:] = field_fourier[i,j,k,0]*xi *TWOPIIMG # vector field from scalar data + grad_fourier[i,j,k,0,:] = field_fourier[i,j,k,0]*xi *TWOPIIMG # vector field from scalar data if dataType == 'vector': - grad_fourier[i,j,k,1,:] = field_fourier[i,j,k,1]*xi *TWOPIIMG # tensor field from vector data + grad_fourier[i,j,k,1,:] = field_fourier[i,j,k,1]*xi *TWOPIIMG # tensor field from vector data grad_fourier[i,j,k,2,:] = field_fourier[i,j,k,2]*xi *TWOPIIMG return np.fft.fftpack.irfftn(grad_fourier,axes=(0,1,2)).reshape([N,3*n]) @@ -79,7 +79,7 @@ parser.set_defaults(coords = 'ipinitialcoord', if options.vector is None and options.scalar is None: parser.error('no data column specified.') -# --- loop over input files ------------------------------------------------------------------------- +# --- loop over input files ------------------------------------------------------------------------ if filenames == []: filenames = [None] @@ -124,7 +124,7 @@ for name in filenames: table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) for type, data in items.iteritems(): for label in data['active']: - table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in xrange(3 * data['dim'])]) # extend ASCII header with new labels # grad increases the field dimension by one + table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in xrange(3 * data['dim'])]) # extend ASCII header with new labels table.head_write() # --------------- figure out size and grid --------------------------------------------------------- @@ -143,7 +143,8 @@ for name in filenames: stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(gradFFT(size[::-1], # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + stack.append(gradFFT(size[::-1], table.data[:,data['column'][i]:data['column'][i]+data['dim']]. reshape([grid[2],grid[1],grid[0]]+data['shape']))) From 36c88cde02a62a0e13961dcf0227907ee93e0a2f Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 17 Mar 2016 04:13:35 +0100 Subject: [PATCH 14/24] updated version information after successful test of v2.0.0-36-g2b7524e --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 27bd816bb..09b18f22c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-31-g6197e83 +v2.0.0-36-g2b7524e From d42f8d5fd4958bbc1015f4db495995e814fbf823 Mon Sep 17 00:00:00 2001 From: Zhuowen Zhao Date: Thu, 17 Mar 2016 17:40:40 -0400 Subject: [PATCH 15/24] added possibility to read real numbers as microstructure indices (not only integers) --- lib/damask/asciitable.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/lib/damask/asciitable.py b/lib/damask/asciitable.py index 4abe4297d..5212ef00c 100644 --- a/lib/damask/asciitable.py +++ b/lib/damask/asciitable.py @@ -517,21 +517,25 @@ class ASCIItable(): # ------------------------------------------------------------------ def microstructure_read(self, - grid): + grid, + type = 'i'): """read microstructure data (from .geom format)""" - N = grid.prod() # expected number of microstructure indices in data - microstructure = np.zeros(N,'i') # initialize as flat array + def datatype(item): + return int(item) if type.lower() == 'i' else float(item) + + N = grid.prod() # expected number of microstructure indices in data + microstructure = np.zeros(N,type) # initialize as flat array i = 0 while i < N and self.data_read(): items = self.data if len(items) > 2: - if items[1].lower() == 'of': items = [int(items[2])]*int(items[0]) - elif items[1].lower() == 'to': items = range(int(items[0]),1+int(items[2])) - else: items = map(int,items) - else: items = map(int,items) + if items[1].lower() == 'of': items = np.ones(datatype(items[0]))*datatype(items[2]) + elif items[1].lower() == 'to': items = np.arange(datatype(items[0]),1+datatype(items[2])) + else: items = map(datatype,items) + else: items = map(datatype,items) - s = min(len(items), N-i) # prevent overflow of microstructure array + s = min(len(items), N-i) # prevent overflow of microstructure array microstructure[i:i+s] = items[:s] i += s From 0f78f1c74713b389462e2aa745a158c2cdcbfd45 Mon Sep 17 00:00:00 2001 From: Aritra Chakraborty Date: Thu, 17 Mar 2016 18:09:22 -0400 Subject: [PATCH 16/24] possibility to treat floating point numbers as microstructure indices --- processing/pre/geom_canvas.py | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/processing/pre/geom_canvas.py b/processing/pre/geom_canvas.py index 882f32c88..ec171f522 100755 --- a/processing/pre/geom_canvas.py +++ b/processing/pre/geom_canvas.py @@ -28,24 +28,32 @@ parser.add_option('-o', '--offset', help = 'a,b,c offset from old to new origin of grid [%default]') parser.add_option('-f', '--fill', dest = 'fill', - type = 'int', metavar = 'int', + type = 'float', metavar = 'float', help = '(background) canvas grain index. "0" selects maximum microstructure index + 1 [%default]') +parser.add_option('--float', + dest = 'real', + action = 'store_true', + help = 'input data is float [%default]') parser.set_defaults(grid = ['0','0','0'], offset = (0,0,0), fill = 0, + real = False, ) (options, filenames) = parser.parse_args() +datatype = 'f' if options.real else 'i' + + # --- loop over input files ------------------------------------------------------------------------- if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name, - buffered = False, labeled = False) + try: table = damask.ASCIItable(name = name, + buffered = False, + labeled = False) except: continue damask.util.report(scriptName,name) @@ -71,7 +79,7 @@ for name in filenames: # --- read data ------------------------------------------------------------------------------------ - microstructure = table.microstructure_read(info['grid']).reshape(info['grid'],order='F') # read microstructure + microstructure = table.microstructure_read(info['grid'],datatype).reshape(info['grid'],order='F') # read microstructure # --- do work ------------------------------------------------------------------------------------ @@ -85,8 +93,8 @@ for name in filenames: else int(n) for o,n in zip(info['grid'],options.grid)],'i') newInfo['grid'] = np.where(newInfo['grid'] > 0, newInfo['grid'],info['grid']) - microstructure_cropped = np.zeros(newInfo['grid'],'i') - microstructure_cropped.fill(options.fill if options.fill > 0 else microstructure.max()+1) + microstructure_cropped = np.zeros(newInfo['grid'],datatype) + microstructure_cropped.fill(options.fill if options.real or options.fill > 0 else microstructure.max()+1) xindex = list(set(xrange(options.offset[0],options.offset[0]+newInfo['grid'][0])) & \ set(xrange(info['grid'][0]))) yindex = list(set(xrange(options.offset[1],options.offset[1]+newInfo['grid'][1])) & \ @@ -152,9 +160,9 @@ for name in filenames: # --- write microstructure information ------------------------------------------------------------ - formatwidth = int(math.floor(math.log10(microstructure_cropped.max())+1)) + format = '%g' if options.real else '%{}i'.format(int(math.floor(math.log10(microstructure_cropped.max())+1))) table.data = microstructure_cropped.reshape((newInfo['grid'][0],newInfo['grid'][1]*newInfo['grid'][2]),order='F').transpose() - table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ') + table.data_writeArray(format,delimiter=' ') # --- output finalization -------------------------------------------------------------------------- From e39441f36bbff3d2dc4d63e42b2492bb6081c16a Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 17 Mar 2016 19:02:30 -0400 Subject: [PATCH 17/24] corrected usage message to tell about string comparison in filter condition --- processing/post/filterTable.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/filterTable.py b/processing/post/filterTable.py index d375785dc..6723d9faf 100755 --- a/processing/post/filterTable.py +++ b/processing/post/filterTable.py @@ -35,7 +35,7 @@ Filter rows according to condition and columns by either white or black listing. Examples: Every odd row if x coordinate is positive -- " #ip.x# >= 0.0 and #_row_#%2 == 1 ). -All rows where label 'foo' equals 'bar' -- " #foo# == \"bar\" " +All rows where label 'foo' equals 'bar' -- " #s#foo# == 'bar' " """, version = scriptID) From 2506667326235ed45fc461160fca211a80e4ef2c Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 18 Mar 2016 04:13:24 +0100 Subject: [PATCH 18/24] updated version information after successful test of v2.0.0-43-ge39441f --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 09b18f22c..690c6cb1d 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-36-g2b7524e +v2.0.0-43-ge39441f From 6865987ea743fc5a2c0f92c42e37766b1ef90dea Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sun, 20 Mar 2016 15:01:47 -0400 Subject: [PATCH 19/24] restored option to output microstructure indices as one-dimensional list (not clear why that got lost!) --- processing/pre/geom_unpack.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/processing/pre/geom_unpack.py b/processing/pre/geom_unpack.py index a07166cdd..85c60d6e5 100755 --- a/processing/pre/geom_unpack.py +++ b/processing/pre/geom_unpack.py @@ -18,6 +18,14 @@ Unpack geometry files containing ranges "a to b" and/or "n of x" multiples (excl """, version = scriptID) +parser.add_option('-1', '--onedimensional', + dest = 'oneD', + action = 'store_true', + help = 'output geom file with one-dimensional data arrangement') + +parser.set_defaults(oneD = False, + ) + (options, filenames) = parser.parse_args() # --- loop over input files ------------------------------------------------------------------------- @@ -69,7 +77,8 @@ for name in filenames: microstructure = table.microstructure_read(info['grid']) # read microstructure formatwidth = int(math.floor(math.log10(microstructure.max())+1)) # efficient number printing format - table.data = microstructure.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose() + table.data = microstructure if options.oneD else \ + microstructure.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose() table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ') #--- output finalization -------------------------------------------------------------------------- From a41cd3df4068af0cd752d72cc060f4d91204f470 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sun, 20 Mar 2016 16:46:35 -0400 Subject: [PATCH 20/24] added "strict" option to reading microstructure from geom file. will now complain if data mismatches with grid. --- lib/damask/asciitable.py | 8 +++++--- processing/pre/geom_check.py | 27 +++++++++++++++------------ 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/lib/damask/asciitable.py b/lib/damask/asciitable.py index 5212ef00c..62d38e4dd 100644 --- a/lib/damask/asciitable.py +++ b/lib/damask/asciitable.py @@ -518,7 +518,8 @@ class ASCIItable(): # ------------------------------------------------------------------ def microstructure_read(self, grid, - type = 'i'): + type = 'i', + strict = False): """read microstructure data (from .geom format)""" def datatype(item): return int(item) if type.lower() == 'i' else float(item) @@ -537,6 +538,7 @@ class ASCIItable(): s = min(len(items), N-i) # prevent overflow of microstructure array microstructure[i:i+s] = items[:s] - i += s + i += len(items) - return microstructure + return (microstructure, i == N and not self.data_read() if strict # check for proper point count and end of file + else microstructure) diff --git a/processing/pre/geom_check.py b/processing/pre/geom_check.py index d607a1f73..518e8c3b0 100755 --- a/processing/pre/geom_check.py +++ b/processing/pre/geom_check.py @@ -54,12 +54,25 @@ for name in filenames: errors = [] if np.any(info['grid'] < 1): errors.append('invalid grid a b c.') if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.') + +#--- read microstructure information -------------------------------------------------------------- + + if options.data: + microstructure,ok = table.microstructure_read(info['grid'],strict = True) # read microstructure + + if ok: + structure = vtk.vtkIntArray() + structure.SetName('Microstructures') + for idx in microstructure: structure.InsertNextValue(idx) + + else: errors.append('mismatch between data and grid dimension.') + if errors != []: damask.util.croak(errors) table.close(dismiss = True) continue -# --- generate VTK rectilinear grid -------------------------------------------------------------------------------- +# --- generate VTK rectilinear grid --------------------------------------------------------------- grid = vtk.vtkRectilinearGrid() grid.SetDimensions([x+1 for x in info['grid']]) @@ -72,18 +85,8 @@ for name in filenames: elif i == 1: grid.SetYCoordinates(temp) elif i == 2: grid.SetZCoordinates(temp) -#--- read microstructure information -------------------------------------------------------------- - if options.data: - microstructure = table.microstructure_read(info['grid']) # read microstructure - - structure = vtk.vtkIntArray() - structure.SetName('Microstructures') - - for idx in microstructure: - structure.InsertNextValue(idx) - - grid.GetCellData().AddArray(structure) + if options.data: grid.GetCellData().AddArray(structure) # --- write data ----------------------------------------------------------------------------------- if name: From 28259b2c46b805b66dcb14d95765feda18e1a7f4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Mar 2016 23:20:58 +0100 Subject: [PATCH 21/24] switched dNeq and dEq to relative tolerance, removed single precision (makes things complicated and was never used anyway) --- code/crystallite.f90 | 46 ++++++++++++++++---------------------------- code/math.f90 | 31 +++++++---------------------- code/prec.f90 | 33 +++++++++++-------------------- 3 files changed, 35 insertions(+), 75 deletions(-) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 53378a70d..71a5e4743 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3569,11 +3569,7 @@ logical function crystallite_integrateStress(& maxticks external :: & -#if(FLOAT==8) dgesv -#elif(FLOAT==4) - sgesv -#endif !* be pessimistic crystallite_integrateStress = .false. @@ -3756,11 +3752,7 @@ logical function crystallite_integrateStress(& - math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333)) dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine work = math_plain33to9(residuumLp) -#if(FLOAT==8) call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp -#elif(FLOAT==4) - call sgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp -#endif if (ierr /= 0_pInt) then #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then @@ -3849,31 +3841,27 @@ logical function crystallite_integrateStress(& math_mul3333xx3333(dT_dFi3333, dFi_dLi3333))) & - math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333)) work = math_plain33to9(residuumLi) -#if(FLOAT==8) call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li -#elif(FLOAT==4) - call sgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li -#endif - if (ierr /= 0_pInt) then + if (ierr /= 0_pInt) then #ifndef _OPENMP - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then - write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', & - el,mesh_element(1,el),ip,ipc - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then - write(6,*) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333)) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',math_transpose33(Li_constitutive) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',math_transpose33(Liguess) - endif + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then + write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', & + el,mesh_element(1,el),ip,ipc + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then + write(6,*) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333)) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',math_transpose33(Li_constitutive) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',math_transpose33(Liguess) endif -#endif - return endif +#endif + return + endif deltaLi = - math_plain9to33(work) endif diff --git a/code/math.f90 b/code/math.f90 index 2f76c1cb5..bf7460062 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -186,10 +186,6 @@ module math halton_seed_set, & i_to_halton, & prime - external :: & - dsyev, & - dgetrf, & - dgetri contains @@ -811,15 +807,13 @@ function math_invSym3333(A) integer(pInt), dimension(6) :: ipiv6 real(pReal), dimension(6,6) :: temp66_Real real(pReal), dimension(6) :: work6 + external :: & + dgetrf, & + dgetri temp66_real = math_Mandel3333to66(A) -#if(FLOAT==8) call dgetrf(6,6,temp66_real,6,ipiv6,ierr) call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) -#elif(FLOAT==4) - call sgetrf(6,6,temp66_real,6,ipiv6,ierr) - call sgetri(6,temp66_real,6,ipiv6,work6,6,ierr) -#endif if (ierr == 0_pInt) then math_invSym3333 = math_Mandel66to3333(temp66_real) else @@ -847,13 +841,8 @@ subroutine math_invert(myDim,A, InvA, error) logical, intent(out) :: error invA = A -#if(FLOAT==8) call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr) call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) -#elif(FLOAT==4) - call sgetrf(myDim,myDim,invA,myDim,ipiv,ierr) - call sgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) -#endif error = merge(.true.,.false., ierr /= 0_pInt) ! http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html end subroutine math_invert @@ -1937,16 +1926,13 @@ subroutine math_eigenValuesVectorsSym(m,values,vectors,error) real(pReal), dimension(size(m,1)), intent(out) :: values real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors logical, intent(out) :: error - integer(pInt) :: info real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f + external :: & + dsyev vectors = m ! copy matrix to input (doubles as output) array -#if(FLOAT==8) call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) -#elif(FLOAT==4) - call ssyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) -#endif error = (info == 0_pInt) end subroutine math_eigenValuesVectorsSym @@ -2135,16 +2121,13 @@ function math_eigenvaluesSym(m) real(pReal), dimension(:,:), intent(in) :: m real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym real(pReal), dimension(size(m,1),size(m,1)) :: vectors - integer(pInt) :: info real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f + external :: & + dsyev vectors = m ! copy matrix to input (doubles as output) array -#if(FLOAT==8) call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) -#elif(FLOAT==4) - call ssyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) -#endif if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN end function math_eigenvaluesSym diff --git a/code/prec.f90 b/code/prec.f90 index 129681fd5..efa0c76a3 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -4,9 +4,9 @@ !> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH -!> @brief setting precision for real and int type depending on makros "FLOAT" and "INT" +!> @brief setting precision for real and int type !> @details setting precision for real and int type and for DAMASK_NaN. Definition is made -!! depending on makros "FLOAT" and "INT" defined during compilation +!! depending on makro "INT" defined during compilation !! for details on NaN see https://software.intel.com/en-us/forums/topic/294680 !-------------------------------------------------------------------------------------------------- module prec @@ -18,18 +18,7 @@ module prec implicit none private -#if (FLOAT==4) -#if defined(Spectral) || defined(FEM) - SPECTRAL SOLVER AND OWN FEM DO NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION -#endif - integer, parameter, public :: pReal = 4 !< floating point single precition (was selected_real_kind(6,37), number with 6 significant digits, up to 1e+-37) -#ifdef __INTEL_COMPILER - real(pReal), parameter, public :: DAMASK_NaN = Z'7F800001' !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran) -#endif -#ifdef __GFORTRAN__ - real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran) -#endif -#elif (FLOAT==8) +#if (FLOAT==8) integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300) #ifdef __INTEL_COMPILER real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran) @@ -188,29 +177,29 @@ logical elemental pure function prec_isNaN(a) end function prec_isNaN - !-------------------------------------------------------------------------------------------------- !> @brief equality comparison for double precision -! replaces "==" but for certain (absolute) tolerance. Counterpart to dNeq +! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq +! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm !-------------------------------------------------------------------------------------------------- logical elemental pure function dEq(a,b,tol) real(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 1.0e-15_pReal - dEq = merge(.True., .False.,abs(a-b) < merge(tol,eps,present(tol))) + real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C + dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function dEq - !-------------------------------------------------------------------------------------------------- !> @brief inequality comparison for double precision -! replaces "!=" but for certain (absolute) tolerance. Counterpart to dEq +! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq +! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm !-------------------------------------------------------------------------------------------------- logical elemental pure function dNeq(a,b,tol) real(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 1.0e-15_pReal - dNeq = merge(.True., .False.,abs(a-b) < merge(tol,eps,present(tol))) + real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C + dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function dNeq end module prec From 01a5343ca20ec2c38449afc8137413b6a7135a62 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sun, 20 Mar 2016 20:58:10 -0400 Subject: [PATCH 22/24] added another backgroundMessage animation. generally cleaned code a slight bit. --- lib/damask/util.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/lib/damask/util.py b/lib/damask/util.py index 9ce42759a..c64a8cd00 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -106,6 +106,7 @@ class backgroundMessage(threading.Thread): choices = {'bounce': ['_', 'o', 'O', u'\u00B0', u'\u203e',u'\u203e',u'\u00B0','O','o','_'], + 'spin': [u'\u25dc',u'\u25dd',u'\u25de',u'\u25df'], 'circle': [u'\u25f4',u'\u25f5',u'\u25f6',u'\u25f7'], 'hexagon': [u'\u2b22',u'\u2b23'], 'square': [u'\u2596',u'\u2598',u'\u259d',u'\u2597'], @@ -229,7 +230,7 @@ def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0, return shape(res), dt def _int2extGrad(p_int, bounds): - """Calculate the gradients of transforming the internal (unconstrained) to external (constained) parameter.""" + """Calculate the gradients of transforming the internal (unconstrained) to external (constrained) parameter.""" grad = np.empty_like(p_int) for i, (x, bound) in enumerate(zip(p_int, bounds)): lower, upper = bound @@ -431,22 +432,20 @@ def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw): else: pcov = np.inf - if return_full: - return popt, pcov, infodict, errmsg, ier - else: - return popt, pcov + return popt, pcov, infodict, errmsg, ier if return_full else popt, pcov -def execute(cmd,streamIn=None,wd='./'): +def execute(cmd,streamIn = None,wd = './'): """executes a command in given directory and returns stdout and stderr for optional stdin""" - initialPath=os.getcwd() + initialPath = os.getcwd() os.chdir(wd) - process = subprocess.Popen(shlex.split(cmd),stdout=subprocess.PIPE,stderr = subprocess.PIPE,stdin=subprocess.PIPE) - if streamIn is not None: - out,error = [i.replace("\x08","") for i in process.communicate(streamIn.read())] - else: - out,error =[i.replace("\x08","") for i in process.communicate()] + process = subprocess.Popen(shlex.split(cmd), + stdout = subprocess.PIPE, + stderr = subprocess.PIPE, + stdin = subprocess.PIPE) + out,error = [i.replace("\x08","") for i in (process.communicate if streamIn is None + else process.communicate(streamIn.read()))] os.chdir(initialPath) - if process.returncode !=0: raise RuntimeError(cmd+' failed with returncode '+str(process.returncode)) + if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) return out,error From 8cd1f027b303b9ca7ff8c8e03bcbcbdc04b2e167 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 21 Mar 2016 08:51:56 -0400 Subject: [PATCH 23/24] bug fixing. --- lib/damask/asciitable.py | 3 +-- lib/damask/util.py | 38 ++++++++++++++++++++------------------ 2 files changed, 21 insertions(+), 20 deletions(-) diff --git a/lib/damask/asciitable.py b/lib/damask/asciitable.py index 62d38e4dd..f76cf6627 100644 --- a/lib/damask/asciitable.py +++ b/lib/damask/asciitable.py @@ -540,5 +540,4 @@ class ASCIItable(): microstructure[i:i+s] = items[:s] i += len(items) - return (microstructure, i == N and not self.data_read() if strict # check for proper point count and end of file - else microstructure) + return (microstructure, i == N and not self.data_read()) if strict else microstructure # check for proper point count and end of file diff --git a/lib/damask/util.py b/lib/damask/util.py index c64a8cd00..244aaa25c 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -58,6 +58,24 @@ def emph(what): """emphasizes string on screen""" return bcolors.BOLD+srepr(what)+bcolors.ENDC +# ----------------------------- +def execute(cmd, + streamIn = None, + wd = './'): + """executes a command in given directory and returns stdout and stderr for optional stdin""" + initialPath = os.getcwd() + os.chdir(wd) + process = subprocess.Popen(shlex.split(cmd), + stdout = subprocess.PIPE, + stderr = subprocess.PIPE, + stdin = subprocess.PIPE) + out,error = [i.replace("\x08","") for i in (process.communicate() if streamIn is None + else process.communicate(streamIn.read()))] + os.chdir(initialPath) + if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) + return out,error + + # ----------------------------- # Matlab like trigonometric functions that take and return angles in degrees. # ----------------------------- @@ -75,7 +93,7 @@ def gridLocation(idx,res): # ----------------------------- def gridIndex(location,res): - return ( location[0] % res[0] + \ + return ( location[0] % res[0] + \ ( location[1] % res[1]) * res[0] + \ ( location[2] % res[2]) * res[1] * res[0] ) @@ -432,20 +450,4 @@ def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw): else: pcov = np.inf - return popt, pcov, infodict, errmsg, ier if return_full else popt, pcov - - -def execute(cmd,streamIn = None,wd = './'): - """executes a command in given directory and returns stdout and stderr for optional stdin""" - initialPath = os.getcwd() - os.chdir(wd) - process = subprocess.Popen(shlex.split(cmd), - stdout = subprocess.PIPE, - stderr = subprocess.PIPE, - stdin = subprocess.PIPE) - out,error = [i.replace("\x08","") for i in (process.communicate if streamIn is None - else process.communicate(streamIn.read()))] - os.chdir(initialPath) - if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) - return out,error - + return (popt, pcov, infodict, errmsg, ier) if return_full else (popt, pcov) From 1b22b6e4ee958d0eb67d4fc3f8e7f2416d184eae Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 21 Mar 2016 10:18:58 -0400 Subject: [PATCH 24/24] replaced variable interpolation into strings by "%" with .format() method --- lib/damask/test/test.py | 145 ++++++++++++++++++++-------------------- 1 file changed, 74 insertions(+), 71 deletions(-) diff --git a/lib/damask/test/test.py b/lib/damask/test/test.py index 523f25d60..826b190c1 100644 --- a/lib/damask/test/test.py +++ b/lib/damask/test/test.py @@ -22,11 +22,11 @@ class Test(): logger = logging.getLogger() logger.setLevel(0) - fh = logging.FileHandler('test.log') # create file handler which logs even debug messages + fh = logging.FileHandler('test.log') # create file handler which logs even debug messages fh.setLevel(logging.DEBUG) full = logging.Formatter('%(asctime)s - %(levelname)s: \n%(message)s') fh.setFormatter(full) - ch = logging.StreamHandler(stream=sys.stdout) # create console handler with a higher log level + ch = logging.StreamHandler(stream=sys.stdout) # create console handler with a higher log level ch.setLevel(logging.INFO) # create formatter and add it to the handlers plain = logging.Formatter('%(message)s') @@ -60,10 +60,10 @@ class Test(): try: self.postprocess(variant) if not self.compare(variant): - return variant+1 # return culprit + return variant+1 # return culprit except Exception as e : - logging.critical('\nWARNING:\n %s\n'%e) - return variant+1 # return culprit + logging.critical('\nWARNING:\n {}\n'.format(e)) + return variant+1 # return culprit return 0 else: if not self.testPossible(): return -1 @@ -74,13 +74,13 @@ class Test(): self.prepare(variant) self.run(variant) self.postprocess(variant) - if self.updateRequested: # update requested + if self.updateRequested: # update requested self.update(variant) - elif not (self.options.accept or self.compare(variant)): # no update, do comparison - return variant+1 # return culprit + elif not (self.options.accept or self.compare(variant)): # no update, do comparison + return variant+1 # return culprit except Exception as e : - logging.critical('\nWARNING:\n %s\n'%e) - return variant+1 # return culprit + logging.critical('\nWARNING:\n {}\n'.format(e)) + return variant+1 # return culprit return 0 def testPossible(self): @@ -94,13 +94,13 @@ class Test(): try: shutil.rmtree(self.dirCurrent()) except: - logging.warning('removal of directory "%s" not possible...'%(self.dirCurrent())) + logging.warning('removal of directory "{}" not possible...'.format(self.dirCurrent())) status = status and False try: os.mkdir(self.dirCurrent()) except: - logging.critical('creation of directory "%s" failed...'%(self.dirCurrent())) + logging.critical('creation of directory "{}" failed...'.format(self.dirCurrent())) status = status and False return status @@ -193,19 +193,19 @@ class Test(): try: shutil.copy2(self.fileInReference(file),self.fileInCurrent(targetfiles[i])) except: - logging.critical('Reference2Current: Unable to copy file %s'%file) + logging.critical('Reference2Current: Unable to copy file "{}"'.format(file)) def copy_Base2Current(self,sourceDir,sourcefiles=[],targetfiles=[]): - source=os.path.normpath(os.path.join(self.dirBase,'../../../'+sourceDir)) + source=os.path.normpath(os.path.join(self.dirBase,'../../..',sourceDir)) if len(targetfiles) == 0: targetfiles = sourcefiles for i,file in enumerate(sourcefiles): try: shutil.copy2(os.path.join(source,file),self.fileInCurrent(targetfiles[i])) except: logging.error(os.path.join(source,file)) - logging.critical('Base2Current: Unable to copy file %s'%file) + logging.critical('Base2Current: Unable to copy file "{}"'.format(file)) def copy_Current2Reference(self,sourcefiles=[],targetfiles=[]): @@ -215,7 +215,7 @@ class Test(): try: shutil.copy2(self.fileInCurrent(file),self.fileInReference(targetfiles[i])) except: - logging.critical('Current2Reference: Unable to copy file %s'%file) + logging.critical('Current2Reference: Unable to copy file "{}"'.format(file)) def copy_Proof2Current(self,sourcefiles=[],targetfiles=[]): @@ -225,7 +225,7 @@ class Test(): try: shutil.copy2(self.fileInProof(file),self.fileInCurrent(targetfiles[i])) except: - logging.critical('Proof2Current: Unable to copy file %s'%file) + logging.critical('Proof2Current: Unable to copy file "{}"'.format(file)) def copy_Current2Current(self,sourcefiles=[],targetfiles=[]): @@ -234,7 +234,7 @@ class Test(): try: shutil.copy2(self.fileInReference(file),self.fileInCurrent(targetfiles[i])) except: - logging.critical('Current2Current: Unable to copy file %s'%file) + logging.critical('Current2Current: Unable to copy file "{}"'.format(file)) def execute_inCurrentDir(self,cmd,streamIn=None): @@ -252,7 +252,7 @@ class Test(): def compare_Array(self,File1,File2): import numpy as np - logging.info('comparing\n '+File1+'\n '+File2) + logging.info('\n '.join(['comparing',File1,File2])) table1 = damask.ASCIItable(name=File1,readonly=True) table1.head_read() len1=len(table1.info)+2 @@ -270,8 +270,9 @@ class Test(): max_loc=np.argmax(abs(refArrayNonZero[curArray.nonzero()]/curArray[curArray.nonzero()]-1.)) refArrayNonZero = refArrayNonZero[curArray.nonzero()] curArray = curArray[curArray.nonzero()] - print(' ********\n * maximum relative error %e for %e and %e\n ********' - %(max_err, refArrayNonZero[max_loc],curArray[max_loc])) + print(' ********\n * maximum relative error {} between {} and {}\n ********'.format(max_err, + refArrayNonZero[max_loc], + curArray[max_loc])) return max_err else: raise Exception('mismatch in array size to compare') @@ -295,7 +296,7 @@ class Test(): absoluteTolerance=False,perLine=False,skipLines=[]): import numpy as np - logging.info('comparing ASCII Tables\n %s \n %s'%(file0,file1)) + logging.info('\n '.join(['comparing ASCII Tables',file0,file1])) if normHeadings == '': normHeadings = headings0 # check if comparison is possible and determine lenght of columns @@ -315,7 +316,7 @@ class Test(): for i in xrange(dataLength): if headings0[i]['shape'] != headings1[i]['shape']: - raise Exception('shape mismatch when comparing %s with %s '%(headings0[i]['label'],headings1[i]['label'])) + raise Exception('shape mismatch between {} and {} '.format(headings0[i]['label'],headings1[i]['label'])) shape[i] = headings0[i]['shape'] for j in xrange(np.shape(shape[i])[0]): length[i] *= shape[i][j] @@ -323,7 +324,9 @@ class Test(): for j in xrange(np.shape(normShape[i])[0]): normLength[i] *= normShape[i][j] else: - raise Exception('trying to compare %i with %i normed by %i data sets'%(len(headings0),len(headings1),len(normHeadings))) + raise Exception('trying to compare {} with {} normed by {} data sets'.format(len(headings0), + len(headings1), + len(normHeadings))) table0 = damask.ASCIItable(name=file0,readonly=True) table0.head_read() @@ -331,37 +334,34 @@ class Test(): table1.head_read() for i in xrange(dataLength): - key0 = {True :'1_%s', - False:'%s' }[length[i]>1]%headings0[i]['label'] - key1 = {True :'1_%s', - False:'%s' }[length[i]>1]%headings1[i]['label'] - normKey = {True :'1_%s', - False:'%s' }[normLength[i]>1]%normHeadings[i]['label'] + key0 = ('1_' if length[i]>1 else '') + headings0[i]['label'] + key1 = ('1_' if length[i]>1 else '') + headings1[i]['label'] + normKey = ('1_' if normLength[i]>1 else '') + normHeadings[i]['label'] if key0 not in table0.labels: - raise Exception('column %s not found in 1. table...\n'%key0) + raise Exception('column {} not found in 1. table...\n'.format(key0)) elif key1 not in table1.labels: - raise Exception('column %s not found in 2. table...\n'%key1) + raise Exception('column {} not found in 2. table...\n'.format(key1)) elif normKey not in table0.labels: - raise Exception('column %s not found in 1. table...\n'%normKey) + raise Exception('column {} not found in 1. table...\n'.format(normKey)) else: - column[0][i] = table0.labels.index(key0) # remember columns of requested data - column[1][i] = table1.labels.index(key1) # remember columns of requested data in second column - normColumn[i] = table0.labels.index(normKey) # remember columns of requested data in second column + column[0][i] = table0.labels.index(key0) + column[1][i] = table1.labels.index(key1) + normColumn[i] = table0.labels.index(normKey) line0 = 0 - while table0.data_read(): # read next data line of ASCII table + while table0.data_read(): # read next data line of ASCII table if line0 not in skipLines: for i in xrange(dataLength): myData = np.array(map(float,table0.data[column[0][i]:\ - column[0][i]+length[i]]),'d') + column[0][i]+length[i]]),'d') normData = np.array(map(float,table0.data[normColumn[i]:\ - normColumn[i]+normLength[i]]),'d') + normColumn[i]+normLength[i]]),'d') data[i] = np.append(data[i],np.reshape(myData,shape[i])) if normType == 'pInf': norm[i] = np.append(norm[i],np.max(np.abs(normData))) else: norm[i] = np.append(norm[i],np.linalg.norm(np.reshape(normData,normShape[i]),normType)) - line0 +=1 + line0 += 1 for i in xrange(dataLength): if not perLine: norm[i] = [np.max(norm[i]) for j in xrange(line0-len(skipLines))] @@ -370,12 +370,12 @@ class Test(): norm[i] = [1.0 for j in xrange(line0-len(skipLines))] absTol[i] = True if perLine: - logging.warning('At least one norm of %s in 1. table is 0.0, using absolute tolerance'%headings0[i]['label']) + logging.warning('At least one norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) else: - logging.warning('Maximum norm of %s in 1. table is 0.0, using absolute tolerance'%headings0[i]['label']) + logging.warning('Maximum norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) line1 = 0 - while table1.data_read(): # read next data line of ASCII table + while table1.data_read(): # read next data line of ASCII table if line1 not in skipLines: for i in xrange(dataLength): myData = np.array(map(float,table1.data[column[1][i]:\ @@ -384,21 +384,25 @@ class Test(): norm[i][line1-len(skipLines)]) line1 +=1 - if (line0 != line1): raise Exception('found %s lines in 1. table and %s in 2. table'%(line0,line1)) + if (line0 != line1): raise Exception('found {} lines in 1. table but {} in 2. table'.format(line0,line1)) logging.info(' ********') for i in xrange(dataLength): if absTol[i]: - logging.info(' * maximum absolute error %e for %s and %s'%(maxError[i],headings0[i]['label'],headings1[i]['label'])) + logging.info(' * maximum absolute error {} between {} and {}'.format(maxError[i], + headings0[i]['label'], + headings1[i]['label'])) else: - logging.info(' * maximum relative error %e for %s and %s'%(maxError[i],headings0[i]['label'],headings1[i]['label'])) + logging.info(' * maximum relative error {} between {} and {}'.format(maxError[i], + headings0[i]['label'], + headings1[i]['label'])) logging.info(' ********') return maxError def compare_TablesStatistically(self, - files = [None,None], # list of file names - columns = [None], # list of list of column labels (per file) + files = [None,None], # list of file names + columns = [None], # list of list of column labels (per file) meanTol = 1.0e-4, stdTol = 1.0e-6, preFilter = 1.0e-9): @@ -407,18 +411,18 @@ class Test(): threshold can be used to ignore small values (a negative number disables this feature) """ - if not (isinstance(files, Iterable) and not isinstance(files, str)): # check whether list of files is requested + if not (isinstance(files, Iterable) and not isinstance(files, str)): # check whether list of files is requested files = [str(files)] tables = [damask.ASCIItable(name = filename,readonly = True) for filename in files] for table in tables: table.head_read() - columns += [columns[0]]*(len(files)-len(columns)) # extend to same length as files - columns = columns[:len(files)] # truncate to same length as files + columns += [columns[0]]*(len(files)-len(columns)) # extend to same length as files + columns = columns[:len(files)] # truncate to same length as files for i,column in enumerate(columns): - if column is None: columns[i] = tables[i].labels # if no column is given, read all + if column is None: columns[i] = tables[i].labels # if no column is given, read all logging.info('comparing ASCIItables statistically') for i in xrange(len(columns)): @@ -428,7 +432,7 @@ class Test(): ) logging.info(files[i]+':'+','.join(columns[i])) - if len(files) < 2: return True # single table is always close to itself... + if len(files) < 2: return True # single table is always close to itself... data = [] for table,labels in zip(tables,columns): @@ -443,16 +447,16 @@ class Test(): normedDelta = np.where(normBy>preFilter,delta/normBy,0.0) mean = np.amax(np.abs(np.mean(normedDelta,0))) std = np.amax(np.std(normedDelta,0)) - logging.info('mean: %f'%mean) - logging.info('std: %f'%std) + logging.info('mean: {:f}'.format(mean)) + logging.info('std: {:f}'.format(std)) return (mean0.0, maximum, 1) # do not devide by zero for empty columns + maximum = np.where(maximum >0.0, maximum, 1) # avoid div by zero for empty columns for i in xrange(len(data)): data[i] /= maximum mask = np.zeros_like(table.data,dtype='bool') for table in data: - mask |= np.where(np.abs(table) 1])) + logging.critical(('The test' if len(self.variants) == 1 else 'All {} tests'.format(len(self.variants))) + ' passed') logging.critical('\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n') return 0 if culprit == -1: logging.warning('Warning: Could not start test') return 0 else: - logging.critical(' ********\n * Test %i failed...\n ********'%(culprit)) + logging.critical(' ********\n * Test {} failed...\n ********'.format(culprit)) logging.critical('\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n') return culprit