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

Merge branch 'development' into MiscImprovements

This commit is contained in:
Martin Diehl 2019-12-21 06:53:56 +01:00
parent 6114f5c325 0fdc880e2c
commit f0d3b29b82
36 changed files with 250 additions and 3986 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
.gitlab-ci.yml
View File

@ -499,18 +499,6 @@ GridSolver:
- master - master
- release - release
Processing:
stage: createDocumentation
script:
- cd $DAMASKROOT/processing/pre
- $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py
- cd $DAMASKROOT/processing/post
- rm vtk2ang.py DAD*.py
- $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py
except:
- master
- release
################################################################################################## ##################################################################################################
backupData: backupData:
stage: saveDocumentation stage: saveDocumentation
@ -521,7 +509,6 @@ backupData:
- mv $TESTROOT/performance/time.png $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/ - mv $TESTROOT/performance/time.png $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
- mv $TESTROOT/performance/memory.png $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/ - mv $TESTROOT/performance/memory.png $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
- mv $DAMASKROOT/PRIVATE/documenting/DAMASK_* $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/ - mv $DAMASKROOT/PRIVATE/documenting/DAMASK_* $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
- mv $DAMASKROOT/processing $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
only: only:
- development - development

2
PRIVATE

@ -1 +1 @@
Subproject commit 524e86c117d816e3bd873eed7663e258a6f2e139 Subproject commit 036faecca39b46fd2328597ca858cbb04e37f79a

2
VERSION
View File

@ -1 +1 @@
v2.0.3-1255-ga6da8fdd v2.0.3-1294-g034367fa

2
examples/ConfigFiles/Homogenization_Thermal_Conduction.config
View File

@ -1,3 +1,3 @@
thermal conduction thermal conduction
initialT 300.0 t0 270.0
(output) temperature (output) temperature

2
examples/SpectralMethod/EshelbyInclusion/material.config
View File

@ -6,7 +6,7 @@
mech none # isostrain 1 grain mech none # isostrain 1 grain
thermal adiabatic # thermal strain (stress) induced mass transport thermal adiabatic # thermal strain (stress) induced mass transport
initialT 300.0 t0 330.0
(output) temperature (output) temperature
#-------------------# #-------------------#

1431
processing/misc/yieldSurface.py
View File

File diff suppressed because it is too large Load Diff

1513
processing/misc/yieldSurfaceFast.py
View File

File diff suppressed because it is too large Load Diff

1
processing/pre/geom_fromDREAM3D.py
View File

@ -145,7 +145,6 @@ for name in filenames:
config_header += ['<microstructure>'] config_header += ['<microstructure>']
for i in range(np.nanmax(microstructure)): for i in range(np.nanmax(microstructure)):
config_header += ['[{}{}]'.format(label,i+1), config_header += ['[{}{}]'.format(label,i+1),
'crystallite 1',
'(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(phase[i],i+1), '(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(phase[i],i+1),
] ]

3
processing/pre/geom_fromOsteonGeometry.py
View File

@ -126,15 +126,12 @@ for i in range(3,np.max(microstructure)):
config_header = ['<microstructure>', config_header = ['<microstructure>',
'[canal]', '[canal]',
'crystallite 1',
'(constituent)\tphase 1\ttexture 1\tfraction 1.0', '(constituent)\tphase 1\ttexture 1\tfraction 1.0',
'[interstitial]', '[interstitial]',
'crystallite 1',
'(constituent)\tphase 2\ttexture 2\tfraction 1.0' '(constituent)\tphase 2\ttexture 2\tfraction 1.0'
] ]
for i in range(3,np.max(microstructure)): for i in range(3,np.max(microstructure)):
config_header += ['[Point{}]'.format(i-2), config_header += ['[Point{}]'.format(i-2),
'crystallite 1',
'(constituent)\tphase 3\ttexture {}\tfraction 1.0'.format(i) '(constituent)\tphase 3\ttexture {}\tfraction 1.0'.format(i)
] ]

1
processing/pre/geom_fromVoronoiTessellation.py
View File

@ -290,7 +290,6 @@ for name in filenames:
config_header += ['<microstructure>'] config_header += ['<microstructure>']
for ID in grainIDs: for ID in grainIDs:
config_header += ['[Grain{}]'.format(ID), config_header += ['[Grain{}]'.format(ID),
'crystallite 1',
'(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(options.phase,ID) '(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(options.phase,ID)
] ]

26
processing/pre/hybridIA_linODFsampling.py
View File

@ -19,7 +19,7 @@ def integerFactorization(i):
return j return j
def binAsBins(bin,intervals): def binAsBins(bin,intervals):
"""Explode compound bin into 3D bins list""" """Explode compound bin into 3D bins list."""
bins = [0]*3 bins = [0]*3
bins[0] = (bin//(intervals[1] * intervals[2])) % intervals[0] bins[0] = (bin//(intervals[1] * intervals[2])) % intervals[0]
bins[1] = (bin//intervals[2]) % intervals[1] bins[1] = (bin//intervals[2]) % intervals[1]
@ -27,17 +27,17 @@ def binAsBins(bin,intervals):
return bins return bins
def binsAsBin(bins,intervals): def binsAsBin(bins,intervals):
"""Implode 3D bins into compound bin""" """Implode 3D bins into compound bin."""
return (bins[0]*intervals[1] + bins[1])*intervals[2] + bins[2] return (bins[0]*intervals[1] + bins[1])*intervals[2] + bins[2]
def EulersAsBins(Eulers,intervals,deltas,center): def EulersAsBins(Eulers,intervals,deltas,center):
"""Return list of Eulers translated into 3D bins list""" """Return list of Eulers translated into 3D bins list."""
return [int((euler+(0.5-center)*delta)//delta)%interval \ return [int((euler+(0.5-center)*delta)//delta)%interval \
for euler,delta,interval in zip(Eulers,deltas,intervals) \ for euler,delta,interval in zip(Eulers,deltas,intervals) \
] ]
def binAsEulers(bin,intervals,deltas,center): def binAsEulers(bin,intervals,deltas,center):
"""Compound bin number translated into list of Eulers""" """Compound bin number translated into list of Eulers."""
Eulers = [0.0]*3 Eulers = [0.0]*3
Eulers[2] = (bin%intervals[2] + center)*deltas[2] Eulers[2] = (bin%intervals[2] + center)*deltas[2]
Eulers[1] = (bin//intervals[2]%intervals[1] + center)*deltas[1] Eulers[1] = (bin//intervals[2]%intervals[1] + center)*deltas[1]
@ -45,7 +45,7 @@ def binAsEulers(bin,intervals,deltas,center):
return Eulers return Eulers
def directInvRepetitions(probability,scale): def directInvRepetitions(probability,scale):
"""Calculate number of samples drawn by direct inversion""" """Calculate number of samples drawn by direct inversion."""
nDirectInv = 0 nDirectInv = 0
for bin in range(len(probability)): # loop over bins for bin in range(len(probability)): # loop over bins
nDirectInv += int(round(probability[bin]*scale)) # calc repetition nDirectInv += int(round(probability[bin]*scale)) # calc repetition
@ -56,7 +56,7 @@ def directInvRepetitions(probability,scale):
# ----- efficient algorithm --------- # ----- efficient algorithm ---------
def directInversion (ODF,nSamples): def directInversion (ODF,nSamples):
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)""" """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
nOptSamples = max(ODF['nNonZero'],nSamples) # random subsampling if too little samples requested nOptSamples = max(ODF['nNonZero'],nSamples) # random subsampling if too little samples requested
nInvSamples = 0 nInvSamples = 0
@ -118,7 +118,7 @@ def directInversion (ODF,nSamples):
# ----- trial and error algorithms --------- # ----- trial and error algorithms ---------
def MonteCarloEulers (ODF,nSamples): def MonteCarloEulers (ODF,nSamples):
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)""" """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
countMC = 0 countMC = 0
maxdV_V = max(ODF['dV_V']) maxdV_V = max(ODF['dV_V'])
orientations = np.zeros((nSamples,3),'f') orientations = np.zeros((nSamples,3),'f')
@ -141,7 +141,7 @@ def MonteCarloEulers (ODF,nSamples):
def MonteCarloBins (ODF,nSamples): def MonteCarloBins (ODF,nSamples):
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)""" """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
countMC = 0 countMC = 0
maxdV_V = max(ODF['dV_V']) maxdV_V = max(ODF['dV_V'])
orientations = np.zeros((nSamples,3),'f') orientations = np.zeros((nSamples,3),'f')
@ -163,7 +163,7 @@ def MonteCarloBins (ODF,nSamples):
def TothVanHoutteSTAT (ODF,nSamples): def TothVanHoutteSTAT (ODF,nSamples):
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)""" """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
orientations = np.zeros((nSamples,3),'f') orientations = np.zeros((nSamples,3),'f')
reconstructedODF = np.zeros(ODF['nBins'],'f') reconstructedODF = np.zeros(ODF['nBins'],'f')
unitInc = 1.0/nSamples unitInc = 1.0/nSamples
@ -211,10 +211,6 @@ parser.add_option('-p','--phase',
dest = 'phase', dest = 'phase',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
help = 'phase index to be used [%default]') help = 'phase index to be used [%default]')
parser.add_option('--crystallite',
dest = 'crystallite',
type = 'int', metavar = 'int',
help = 'crystallite index to be used [%default]')
parser.add_option('-r', '--rnd', parser.add_option('-r', '--rnd',
dest = 'randomSeed', dest = 'randomSeed',
type = 'int', metavar = 'int', \ type = 'int', metavar = 'int', \
@ -223,7 +219,6 @@ parser.set_defaults(randomSeed = None,
number = 500, number = 500,
algorithm = 'IA', algorithm = 'IA',
phase = 1, phase = 1,
crystallite = 1,
ang = True, ang = True,
) )
@ -240,7 +235,7 @@ if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try:
table = damask.ASCIItable(name = name, buffered = False, readonly=True) table = damask.ASCIItable(name = name, buffered = False, readonly=True)
except: except IOError:
continue continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
@ -351,7 +346,6 @@ for name in filenames:
for i,ID in enumerate(range(nSamples)): for i,ID in enumerate(range(nSamples)):
materialConfig += ['[Grain%s]'%(str(ID+1).zfill(formatwidth)), materialConfig += ['[Grain%s]'%(str(ID+1).zfill(formatwidth)),
'crystallite %i'%options.crystallite,
'(constituent) phase %i texture %s fraction 1.0'%(options.phase,str(ID+1).rjust(formatwidth)), '(constituent) phase %i texture %s fraction 1.0'%(options.phase,str(ID+1).rjust(formatwidth)),
] ]

16
processing/pre/patchFromReconstructedBoundaries.py
View File

@ -78,13 +78,11 @@ def rcbOrientationParser(content,idcolumn):
damask.util.croak('You might not have chosen the correct column for the grain IDs! '+ damask.util.croak('You might not have chosen the correct column for the grain IDs! '+
'Please check the "--id" option.') 'Please check the "--id" option.')
raise raise
except:
raise
return grains return grains
def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn): def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
"""Parser for TSL-OIM reconstructed boundary files""" """Parser for TSL-OIM reconstructed boundary files."""
# find bounding box # find bounding box
boxX = [1.*sys.maxint,-1.*sys.maxint] boxX = [1.*sys.maxint,-1.*sys.maxint]
boxY = [1.*sys.maxint,-1.*sys.maxint] boxY = [1.*sys.maxint,-1.*sys.maxint]
@ -99,8 +97,6 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
damask.util.croak('You might not have chosen the correct column for the segment end points! '+ damask.util.croak('You might not have chosen the correct column for the segment end points! '+
'Please check the "--segment" option.') 'Please check the "--segment" option.')
raise raise
except:
raise
(x[0],y[0]) = (M[0]*x[0]+M[1]*y[0],M[2]*x[0]+M[3]*y[0]) # apply transformation to coordinates (x[0],y[0]) = (M[0]*x[0]+M[1]*y[0],M[2]*x[0]+M[3]*y[0]) # apply transformation to coordinates
(x[1],y[1]) = (M[0]*x[1]+M[1]*y[1],M[2]*x[1]+M[3]*y[1]) # to get rcb --> Euler system (x[1],y[1]) = (M[0]*x[1]+M[1]*y[1],M[2]*x[1]+M[3]*y[1]) # to get rcb --> Euler system
boxX[0] = min(boxX[0],x[0],x[1]) boxX[0] = min(boxX[0],x[0],x[1])
@ -728,7 +724,7 @@ def image(name,imgsize,marginX,marginY,rcData):
# ------------------------- # -------------------------
def inside(x,y,points): def inside(x,y,points):
"""Tests whether point(x,y) is within polygon described by points""" """Tests whether point(x,y) is within polygon described by points."""
inside = False inside = False
npoints=len(points) npoints=len(points)
(x1,y1) = points[npoints-1] # start with last point of points (x1,y1) = points[npoints-1] # start with last point of points
@ -750,7 +746,7 @@ def inside(x,y,points):
# ------------------------- # -------------------------
def fftbuild(rcData,height,xframe,yframe,grid,extrusion): def fftbuild(rcData,height,xframe,yframe,grid,extrusion):
"""Build array of grain numbers""" """Build array of grain numbers."""
maxX = -1.*sys.maxint maxX = -1.*sys.maxint
maxY = -1.*sys.maxint maxY = -1.*sys.maxint
for line in rcData['point']: # find data range for line in rcData['point']: # find data range
@ -883,7 +879,7 @@ try:
boundaryFile = open(args[0]) boundaryFile = open(args[0])
boundarySegments = boundaryFile.readlines() boundarySegments = boundaryFile.readlines()
boundaryFile.close() boundaryFile.close()
except: except IOError:
damask.util.croak('unable to read boundary file "{}".'.format(args[0])) damask.util.croak('unable to read boundary file "{}".'.format(args[0]))
raise raise
@ -941,19 +937,15 @@ if any(output in options.output for output in ['spectral','mentat']):
for i,grain in enumerate(rcData['grainMapping']): for i,grain in enumerate(rcData['grainMapping']):
config+=['[grain{}]'.format(grain), config+=['[grain{}]'.format(grain),
'crystallite\t1',
'(constituent)\tphase 1\ttexture {}\tfraction 1.0'.format(i+1)] '(constituent)\tphase 1\ttexture {}\tfraction 1.0'.format(i+1)]
if (options.xmargin > 0.0): if (options.xmargin > 0.0):
config+=['[x-margin]', config+=['[x-margin]',
'crystallite\t1',
'(constituent)\tphase 2\ttexture {}\tfraction 1.0\n'.format(len(rcData['grainMapping'])+1)] '(constituent)\tphase 2\ttexture {}\tfraction 1.0\n'.format(len(rcData['grainMapping'])+1)]
if (options.ymargin > 0.0): if (options.ymargin > 0.0):
config+=['[y-margin]', config+=['[y-margin]',
'crystallite\t1',
'(constituent)\tphase 2\ttexture {}\tfraction 1.0\n'.format(len(rcData['grainMapping'])+1)] '(constituent)\tphase 2\ttexture {}\tfraction 1.0\n'.format(len(rcData['grainMapping'])+1)]
if (options.xmargin > 0.0 and options.ymargin > 0.0): if (options.xmargin > 0.0 and options.ymargin > 0.0):
config+=['[xy-margin]', config+=['[xy-margin]',
'crystallite\t1',
'(constituent)\tphase 2\ttexture {}\tfraction 1.0\n'.format(len(rcData['grainMapping'])+1)] '(constituent)\tphase 2\ttexture {}\tfraction 1.0\n'.format(len(rcData['grainMapping'])+1)]
if (options.xmargin > 0.0 or options.ymargin > 0.0): if (options.xmargin > 0.0 or options.ymargin > 0.0):

4
python/damask/dadf5.py
View File

@ -366,7 +366,7 @@ class DADF5():
f[k] f[k]
path.append(k) path.append(k)
except KeyError as e: except KeyError as e:
print('unable to locate geometry dataset: {}'.format(str(e))) pass
for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']):
for oo in self.iter_visible(o): for oo in self.iter_visible(o):
for pp in self.iter_visible(p): for pp in self.iter_visible(p):
@ -375,7 +375,7 @@ class DADF5():
f[k] f[k]
path.append(k) path.append(k)
except KeyError as e: except KeyError as e:
print('unable to locate {} dataset: {}'.format(o,str(e))) pass
return path return path

5
python/damask/table.py
View File

@ -203,7 +203,7 @@ class Table():
'' if info is None else ': {}'.format(info), '' if info is None else ': {}'.format(info),
)) ))
self.shapes[label_new] = self.shapes.pop(label_old) self.shapes = {(label if label is not label_old else label_new):self.shapes[label] for label in self.shapes}
def sort_by(self,labels,ascending=True): def sort_by(self,labels,ascending=True):
@ -234,8 +234,9 @@ class Table():
Filename or file for reading. Filename or file for reading.
""" """
seen = set()
labels = [] labels = []
for l in self.shapes: for l in [x for x in self.data.columns if not (x in seen or seen.add(x))]:
if(self.shapes[l] == (1,)): if(self.shapes[l] == (1,)):
labels.append('{}'.format(l)) labels.append('{}'.format(l))
elif(len(self.shapes[l]) == 1): elif(len(self.shapes[l]) == 1):

257
python/damask/util.py
View File

@ -221,263 +221,6 @@ class return_message():
return srepr(self.message) return srepr(self.message)
def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0,
col_deriv=0, ftol=1.49012e-8, xtol=1.49012e-8,
gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None):
from scipy.optimize import _minpack
"""
Non-linear least square fitting (Levenberg-Marquardt method) with
bounded parameters.
the codes of transformation between int <-> ext refers to the work of
Jonathan J. Helmus: https://github.com/jjhelmus/leastsqbound-scipy
other codes refer to the source code of minpack.py:
An internal parameter list is used to enforce contraints on the fitting
parameters. The transfomation is based on that of MINUIT package.
please see: F. James and M. Winkler. MINUIT User's Guide, 2004.
bounds : list
(min, max) pairs for each parameter, use None for 'min' or 'max'
when there is no bound in that direction.
For example: if there are two parameters needed to be fitting, then
bounds is [(min1,max1), (min2,max2)]
This function is based on 'leastsq' of minpack.py, the annotation of
other parameters can be found in 'least_squares.py'.
"""
def _check_func(checker, argname, thefunc, x0, args, numinputs,
output_shape=None):
from numpy import shape
"""The same as that of minpack.py"""
res = np.atleast_1d(thefunc(*((x0[:numinputs],) + args)))
if (output_shape is not None) and (shape(res) != output_shape):
if (output_shape[0] != 1):
if len(output_shape) > 1:
if output_shape[1] == 1:
return shape(res)
msg = "%s: there is a mismatch between the input and output " \
"shape of the '%s' argument" % (checker, argname)
func_name = getattr(thefunc, '__name__', None)
if func_name:
msg += " '%s'." % func_name
else:
msg += "."
raise TypeError(msg)
if np.issubdtype(res.dtype, np.inexact):
dt = res.dtype
else:
dt = dtype(float)
return shape(res), dt
def _int2extGrad(p_int, bounds):
"""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
if lower is None and upper is None: # No constraints
grad[i] = 1.0
elif upper is None: # only lower bound
grad[i] = x/np.sqrt(x*x + 1.0)
elif lower is None: # only upper bound
grad[i] = -x/np.sqrt(x*x + 1.0)
else: # lower and upper bounds
grad[i] = (upper - lower)*np.cos(x)/2.0
return grad
def _int2extFunc(bounds):
"""Transform internal parameters into external parameters."""
local = [_int2extLocal(b) for b in bounds]
def _transform_i2e(p_int):
p_ext = np.empty_like(p_int)
p_ext[:] = [i(j) for i, j in zip(local, p_int)]
return p_ext
return _transform_i2e
def _ext2intFunc(bounds):
"""Transform external parameters into internal parameters."""
local = [_ext2intLocal(b) for b in bounds]
def _transform_e2i(p_ext):
p_int = np.empty_like(p_ext)
p_int[:] = [i(j) for i, j in zip(local, p_ext)]
return p_int
return _transform_e2i
def _int2extLocal(bound):
"""Transform a single internal parameter to an external parameter."""
lower, upper = bound
if lower is None and upper is None: # no constraints
return lambda x: x
elif upper is None: # only lower bound
return lambda x: lower - 1.0 + np.sqrt(x*x + 1.0)
elif lower is None: # only upper bound
return lambda x: upper + 1.0 - np.sqrt(x*x + 1.0)
else:
return lambda x: lower + ((upper - lower)/2.0)*(np.sin(x) + 1.0)
def _ext2intLocal(bound):
"""Transform a single external parameter to an internal parameter."""
lower, upper = bound
if lower is None and upper is None: # no constraints
return lambda x: x
elif upper is None: # only lower bound
return lambda x: np.sqrt((x - lower + 1.0)**2 - 1.0)
elif lower is None: # only upper bound
return lambda x: np.sqrt((x - upper - 1.0)**2 - 1.0)
else:
return lambda x: np.arcsin((2.0*(x - lower)/(upper - lower)) - 1.0)
i2e = _int2extFunc(bounds)
e2i = _ext2intFunc(bounds)
x0 = np.asarray(x0).flatten()
n = len(x0)
if len(bounds) != n:
raise ValueError('the length of bounds is inconsistent with the number of parameters ')
if not isinstance(args, tuple):
args = (args,)
shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
m = shape[0]
if n > m:
raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
if epsfcn is None:
epsfcn = np.finfo(dtype).eps
def funcWarp(x, *args):
return func(i2e(x), *args)
xi0 = e2i(x0)
if Dfun is None:
if maxfev == 0:
maxfev = 200*(n + 1)
retval = _minpack._lmdif(funcWarp, xi0, args, full_output, ftol, xtol,
gtol, maxfev, epsfcn, factor, diag)
else:
if col_deriv:
_check_func('leastsq', 'Dfun', Dfun, x0, args, n, (n, m))
else:
_check_func('leastsq', 'Dfun', Dfun, x0, args, n, (m, n))
if maxfev == 0:
maxfev = 100*(n + 1)
def DfunWarp(x, *args):
return Dfun(i2e(x), *args)
retval = _minpack._lmder(funcWarp, DfunWarp, xi0, args, full_output, col_deriv,
ftol, xtol, gtol, maxfev, factor, diag)
errors = {0: ["Improper input parameters.", TypeError],
1: ["Both actual and predicted relative reductions "
"in the sum of squares\n are at most %f" % ftol, None],
2: ["The relative error between two consecutive "
"iterates is at most %f" % xtol, None],
3: ["Both actual and predicted relative reductions in "
"the sum of squares\n are at most %f and the "
"relative error between two consecutive "
"iterates is at \n most %f" % (ftol, xtol), None],
4: ["The cosine of the angle between func(x) and any "
"column of the\n Jacobian is at most %f in "
"absolute value" % gtol, None],
5: ["Number of calls to function has reached "
"maxfev = %d." % maxfev, ValueError],
6: ["ftol=%f is too small, no further reduction "
"in the sum of squares\n is possible.""" % ftol,
ValueError],
7: ["xtol=%f is too small, no further improvement in "
"the approximate\n solution is possible." % xtol,
ValueError],
8: ["gtol=%f is too small, func(x) is orthogonal to the "
"columns of\n the Jacobian to machine "
"precision." % gtol, ValueError],
'unknown': ["Unknown error.", TypeError]}
info = retval[-1] # The FORTRAN return value
if info not in [1, 2, 3, 4] and not full_output:
if info in [5, 6, 7, 8]:
np.warnings.warn(errors[info][0], RuntimeWarning)
else:
try:
raise errors[info][1](errors[info][0])
except KeyError:
raise errors['unknown'][1](errors['unknown'][0])
mesg = errors[info][0]
x = i2e(retval[0])
if full_output:
grad = _int2extGrad(retval[0], bounds)
retval[1]['fjac'] = (retval[1]['fjac'].T / np.take(grad,
retval[1]['ipvt'] - 1)).T
cov_x = None
if info in [1, 2, 3, 4]:
from numpy.dual import inv
from numpy.linalg import LinAlgError
perm = np.take(np.eye(n), retval[1]['ipvt'] - 1, 0)
r = np.triu(np.transpose(retval[1]['fjac'])[:n, :])
R = np.dot(r, perm)
try:
cov_x = inv(np.dot(np.transpose(R), R))
except LinAlgError as inverror:
print(inverror)
pass
return (x, cov_x) + retval[1:-1] + (mesg, info)
else:
return (x, info)
def _general_function(params, ydata, xdata, function):
return function(xdata, *params) - ydata
def _weighted_general_function(params, ydata, xdata, function, weights):
return (function(xdata, *params) - ydata)*weights
def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw):
"""Similar as 'curve_fit' in minpack.py."""
if p0 is None:
# determine number of parameters by inspecting the function
import inspect
args, varargs, varkw, defaults = inspect.getargspec(f)
if len(args) < 2:
msg = "Unable to determine number of fit parameters."
raise ValueError(msg)
if 'self' in args:
p0 = [1.0] * (len(args)-2)
else:
p0 = [1.0] * (len(args)-1)
if np.isscalar(p0):
p0 = np.array([p0])
args = (ydata, xdata, f)
if sigma is None:
func = _general_function
else:
func = _weighted_general_function
args += (1.0/np.asarray(sigma),)
return_full = kw.pop('full_output', False)
res = leastsqBound(func, p0, args=args, bounds = bounds, full_output=True, **kw)
(popt, pcov, infodict, errmsg, ier) = res
if ier not in [1, 2, 3, 4]:
msg = "Optimal parameters not found: " + errmsg
raise RuntimeError(msg)
if (len(ydata) > len(p0)) and pcov is not None:
s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq
else:
pcov = np.inf
return (popt, pcov, infodict, errmsg, ier) if return_full else (popt, pcov)
class ThreadPool: class ThreadPool:
"""Pool of threads consuming tasks from a queue.""" """Pool of threads consuming tasks from a queue."""

1
python/setup.py
View File

@ -15,6 +15,7 @@ setuptools.setup(
packages=setuptools.find_packages(), packages=setuptools.find_packages(),
include_package_data=True, include_package_data=True,
install_requires = [ install_requires = [
"pandas",
"scipy", "scipy",
"h5py", "h5py",
"vtk" "vtk"

12
python/tests/test_Table.py
View File

@ -65,11 +65,13 @@ class TestTable:
default.add('nine',d,'random data') default.add('nine',d,'random data')
assert np.allclose(d,default.get('nine')) assert np.allclose(d,default.get('nine'))
def test_rename_equivalent(self,default): def test_rename_equivalent(self):
v = default.get('v') x = np.random.random((5,13))
default.rename('v','u') t = Table(x,{'F':(3,3),'v':(3,),'s':(1,)},['random test data'])
u = default.get('u') s = t.get('s')
assert np.all(v == u) t.rename('s','u')
u = t.get('u')
assert np.all(s == u)
def test_rename_gone(self,default): def test_rename_gone(self,default):
default.rename('v','V') default.rename('v','V')

1
src/CPFEM2.f90
View File

@ -65,7 +65,6 @@ subroutine CPFEM_initAll
call constitutive_init call constitutive_init
call crystallite_init call crystallite_init
call homogenization_init call homogenization_init
call materialpoint_postResults
call CPFEM_init call CPFEM_init
end subroutine CPFEM_initAll end subroutine CPFEM_initAll

5
src/DAMASK_abaqus.f
View File

@ -143,9 +143,6 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
outdatedByNewInc, & outdatedByNewInc, &
outdatedFFN1, & outdatedFFN1, &
lastStep lastStep
use homogenization, only: &
materialpoint_sizeResults, &
materialpoint_results
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -332,7 +329,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
ddsdde(6,:) = ddsdde_h(5,:) ddsdde(6,:) = ddsdde_h(5,:)
end if end if
statev = materialpoint_results(1:min(nstatv,materialpoint_sizeResults),npt,mesh_FEasCP('elem', noel)) statev = 0
if (terminallyIll) pnewdt = 0.5_pReal ! force cutback directly ? if (terminallyIll) pnewdt = 0.5_pReal ! force cutback directly ?
!$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value !$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value

23
src/IO.f90
View File

@ -32,8 +32,7 @@ module IO
IO_intValue, & IO_intValue, &
IO_lc, & IO_lc, &
IO_error, & IO_error, &
IO_warning, & IO_warning
IO_intOut
#if defined(Marc4DAMASK) || defined(Abaqus) #if defined(Marc4DAMASK) || defined(Abaqus)
public :: & public :: &
IO_open_inputFile, & IO_open_inputFile, &
@ -542,26 +541,6 @@ pure function IO_lc(string)
end function IO_lc end function IO_lc
!--------------------------------------------------------------------------------------------------
!> @brief returns format string for integer values without leading zeros
!> @details deprecated, use '(i0)' format specifier
!--------------------------------------------------------------------------------------------------
pure function IO_intOut(intToPrint)
integer, intent(in) :: intToPrint
character(len=41) :: IO_intOut
integer :: N_digits
character(len=19) :: width ! maximum digits for 64 bit integer
character(len=20) :: min_width ! longer for negative values
N_digits = 1 + int(log10(real(max(abs(intToPrint),1))))
write(width, '(I19.19)') N_digits
write(min_width, '(I20.20)') N_digits + merge(1,0,intToPrint < 0)
IO_intOut = 'I'//trim(min_width)//'.'//trim(width)
end function IO_intOut
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx !> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx
!> in ABAQUS either time step is reduced or execution terminated !> in ABAQUS either time step is reduced or execution terminated

110
src/constitutive.f90
View File

@ -37,7 +37,6 @@ module constitutive
integer, public, protected :: & integer, public, protected :: &
constitutive_plasticity_maxSizeDotState, & constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizePostResults, &
constitutive_source_maxSizeDotState constitutive_source_maxSizeDotState
public :: & public :: &
@ -50,7 +49,6 @@ module constitutive
constitutive_SandItsTangents, & constitutive_SandItsTangents, &
constitutive_collectDotState, & constitutive_collectDotState, &
constitutive_collectDeltaState, & constitutive_collectDeltaState, &
constitutive_postResults, &
constitutive_results constitutive_results
contains contains
@ -61,17 +59,9 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_init subroutine constitutive_init
integer, parameter :: FILEUNIT = 204
integer :: & integer :: &
o, & !< counter in output loop
ph, & !< counter in phase loop ph, & !< counter in phase loop
s, & !< counter in source loop s !< counter in source loop
ins !< instance of plasticity/source
integer, dimension(:,:), pointer :: thisSize
character(len=64), dimension(:,:), pointer :: thisOutput
character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready
logical :: knownSource
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialized plasticity ! initialized plasticity
@ -101,58 +91,10 @@ subroutine constitutive_init
if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init
if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init
write(6,'(/,a)') ' <<<+- constitutive init -+>>>' write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6)
mainProcess: if (worldrank == 0) then
!--------------------------------------------------------------------------------------------------
! write description file for constitutive output
call IO_write_jobFile(FILEUNIT,'outputConstitutive')
PhaseLoop: do ph = 1,material_Nphase
activePhase: if (any(material_phaseAt == ph)) then
write(FILEUNIT,'(/,a,/)') '['//trim(config_name_phase(ph))//']'
SourceLoop: do s = 1, phase_Nsources(ph)
knownSource = .true. ! assume valid
sourceType: select case (phase_source(s,ph))
case (SOURCE_damage_isoBrittle_ID) sourceType
ins = source_damage_isoBrittle_instance(ph)
outputName = SOURCE_damage_isoBrittle_label
thisOutput => source_damage_isoBrittle_output
thisSize => source_damage_isoBrittle_sizePostResult
case (SOURCE_damage_isoDuctile_ID) sourceType
ins = source_damage_isoDuctile_instance(ph)
outputName = SOURCE_damage_isoDuctile_label
thisOutput => source_damage_isoDuctile_output
thisSize => source_damage_isoDuctile_sizePostResult
case (SOURCE_damage_anisoBrittle_ID) sourceType
ins = source_damage_anisoBrittle_instance(ph)
outputName = SOURCE_damage_anisoBrittle_label
thisOutput => source_damage_anisoBrittle_output
thisSize => source_damage_anisoBrittle_sizePostResult
case (SOURCE_damage_anisoDuctile_ID) sourceType
ins = source_damage_anisoDuctile_instance(ph)
outputName = SOURCE_damage_anisoDuctile_label
thisOutput => source_damage_anisoDuctile_output
thisSize => source_damage_anisoDuctile_sizePostResult
case default sourceType
knownSource = .false.
end select sourceType
if (knownSource) then
write(FILEUNIT,'(a)') '(source)'//char(9)//trim(outputName)
OutputSourceLoop: do o = 1,size(thisOutput(:,ins))
if(len_trim(thisOutput(o,ins)) > 0) &
write(FILEUNIT,'(a,i4)') trim(thisOutput(o,ins))//char(9),thisSize(o,ins)
enddo OutputSourceLoop
endif
enddo SourceLoop
endif activePhase
enddo PhaseLoop
close(FILEUNIT)
endif mainProcess
constitutive_plasticity_maxSizeDotState = 0 constitutive_plasticity_maxSizeDotState = 0
constitutive_source_maxSizeDotState = 0 constitutive_source_maxSizeDotState = 0
constitutive_source_maxSizePostResults = 0
PhaseLoop2:do ph = 1,material_Nphase PhaseLoop2:do ph = 1,material_Nphase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -169,11 +111,8 @@ subroutine constitutive_init
plasticState(ph)%sizeDotState) plasticState(ph)%sizeDotState)
constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, &
maxval(sourceState(ph)%p(:)%sizeDotState)) maxval(sourceState(ph)%p(:)%sizeDotState))
constitutive_source_maxSizePostResults = max(constitutive_source_maxSizePostResults, &
maxval(sourceState(ph)%p(:)%sizePostResults))
enddo PhaseLoop2 enddo PhaseLoop2
end subroutine constitutive_init end subroutine constitutive_init
@ -639,51 +578,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
end subroutine constitutive_collectDeltaState end subroutine constitutive_collectDeltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results
!--------------------------------------------------------------------------------------------------
function constitutive_postResults(S, Fi, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(sum(sourceState(material_phaseAt(ipc,el))%p(:)%sizePostResults)) :: &
constitutive_postResults
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(in), dimension(3,3) :: &
S !< 2nd Piola Kirchhoff stress
integer :: &
startPos, endPos
integer :: &
i, of, instance !< counter in source loop
constitutive_postResults = 0.0_pReal
endPos = 0
SourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el))
startPos = endPos + 1
endPos = endPos + sourceState(material_phaseAt(ipc,el))%p(i)%sizePostResults
of = material_phasememberAt(ipc,ip,el)
sourceType: select case (phase_source(i,material_phaseAt(ipc,el)))
case (SOURCE_damage_isoBrittle_ID) sourceType
constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(material_phaseAt(ipc,el),of)
case (SOURCE_damage_isoDuctile_ID) sourceType
constitutive_postResults(startPos:endPos) = source_damage_isoDuctile_postResults(material_phaseAt(ipc,el),of)
case (SOURCE_damage_anisoBrittle_ID) sourceType
constitutive_postResults(startPos:endPos) = source_damage_anisoBrittle_postResults(material_phaseAt(ipc,el),of)
case (SOURCE_damage_anisoDuctile_ID) sourceType
constitutive_postResults(startPos:endPos) = source_damage_anisoDuctile_postResults(material_phaseAt(ipc,el),of)
end select sourceType
enddo SourceLoop
end function constitutive_postResults
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes constitutive results to HDF5 output file !> @brief writes constitutive results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

40
src/crystallite.f90
View File

@ -108,7 +108,6 @@ module crystallite
crystallite_stressTangent, & crystallite_stressTangent, &
crystallite_orientations, & crystallite_orientations, &
crystallite_push33ToRef, & crystallite_push33ToRef, &
crystallite_postResults, &
crystallite_results crystallite_results
contains contains
@ -119,7 +118,6 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_init subroutine crystallite_init
integer, parameter :: FILEUNIT=434
logical, dimension(:,:), allocatable :: devNull logical, dimension(:,:), allocatable :: devNull
integer :: & integer :: &
c, & !< counter in integration point component loop c, & !< counter in integration point component loop
@ -233,13 +231,6 @@ subroutine crystallite_init
#endif #endif
enddo enddo
!--------------------------------------------------------------------------------------------------
! write description file for crystallite output
if (worldrank == 0) then
call IO_write_jobFile(FILEUNIT,'outputCrystallite')
write(FILEUNIT,'(/,a,/)') '[not supported anymore]'
close(FILEUNIT)
endif
call config_deallocate('material.config/phase') call config_deallocate('material.config/phase')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -732,37 +723,6 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33)
end function crystallite_push33ToRef end function crystallite_push33ToRef
!--------------------------------------------------------------------------------------------------
!> @brief return results of particular grain
!--------------------------------------------------------------------------------------------------
function crystallite_postResults(ipc, ip, el)
integer, intent(in):: &
el, & !< element index
ip, & !< integration point index
ipc !< grain index
real(pReal), dimension(1+ &
1+sum(sourceState(material_phaseAt(ipc,el))%p(:)%sizePostResults)) :: &
crystallite_postResults
integer :: &
c
crystallite_postResults = 0.0_pReal
crystallite_postResults(1) = 0.0_pReal ! header-like information (length)
c = 1
crystallite_postResults(c+1) = real(sum(sourceState(material_phaseAt(ipc,el))%p(:)%sizePostResults),pReal) ! size of constitutive results
c = c + 1
if (size(crystallite_postResults)-c > 0) &
crystallite_postResults(c+1:size(crystallite_postResults)) = &
constitutive_postResults(crystallite_S(1:3,1:3,ipc,ip,el), crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc, ip, el)
end function crystallite_postResults
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes crystallite results to HDF5 output file !> @brief writes crystallite results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

39
src/damage_local.f90
View File

@ -11,12 +11,11 @@ module damage_local
use source_damage_isoDuctile use source_damage_isoDuctile
use source_damage_anisoBrittle use source_damage_anisoBrittle
use source_damage_anisoDuctile use source_damage_anisoDuctile
use results
implicit none implicit none
private private
integer, dimension(:,:), allocatable, target, public :: &
damage_local_sizePostResult
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
damage_local_output damage_local_output
integer, dimension(:), allocatable, target, public :: & integer, dimension(:), allocatable, target, public :: &
@ -42,7 +41,7 @@ module damage_local
public :: & public :: &
damage_local_init, & damage_local_init, &
damage_local_updateState, & damage_local_updateState, &
damage_local_postResults damage_local_Results
contains contains
@ -66,7 +65,6 @@ subroutine damage_local_init
maxNinstance = count(damage_type == DAMAGE_local_ID) maxNinstance = count(damage_type == DAMAGE_local_ID)
if (maxNinstance == 0) return if (maxNinstance == 0) return
allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_local_output (maxval(homogenization_Noutput),maxNinstance)) allocate(damage_local_output (maxval(homogenization_Noutput),maxNinstance))
damage_local_output = '' damage_local_output = ''
allocate(damage_local_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) allocate(damage_local_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
@ -90,7 +88,6 @@ subroutine damage_local_init
case ('damage') case ('damage')
damage_local_output(i,damage_typeInstance(h)) = outputs(i) damage_local_output(i,damage_typeInstance(h)) = outputs(i)
damage_local_Noutput(instance) = damage_local_Noutput(instance) + 1 damage_local_Noutput(instance) = damage_local_Noutput(instance) + 1
damage_local_sizePostResult(i,damage_typeInstance(h)) = 1
prm%outputID = [prm%outputID , damage_ID] prm%outputID = [prm%outputID , damage_ID]
end select end select
@ -106,7 +103,6 @@ subroutine damage_local_init
! allocate state arrays ! allocate state arrays
sizeState = 1 sizeState = 1
damageState(homog)%sizeState = sizeState damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_local_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog)) allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog)) allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog)) allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog))
@ -211,35 +207,30 @@ end subroutine damage_local_getSourceAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of damage results !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function damage_local_postResults(ip,el) subroutine damage_local_results(homog,group)
integer, intent(in) :: & integer, intent(in) :: homog
ip, & !< integration point character(len=*), intent(in) :: group
el !< element #if defined(PETSc) || defined(DAMASK_HDF5)
real(pReal), dimension(sum(damage_local_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: & integer :: o, instance
damage_local_postResults
integer :: instance, homog, offset, o, c
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog) instance = damage_typeInstance(homog)
associate(prm => param(instance)) associate(prm => param(instance))
c = 0
outputsLoop: do o = 1,size(prm%outputID) outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o)) select case(prm%outputID(o))
case (damage_ID) case (damage_ID)
damage_local_postResults(c+1) = damage(homog)%p(offset) call results_writeDataset(group,damage(homog)%p,'phi',&
c = c + 1 'damage indicator','-')
end select end select
enddo outputsLoop enddo outputsLoop
end associate end associate
#endif
end subroutine damage_local_results
end function damage_local_postResults
end module damage_local end module damage_local

42
src/damage_nonlocal.f90
View File

@ -14,12 +14,11 @@ module damage_nonlocal
use source_damage_isoDuctile use source_damage_isoDuctile
use source_damage_anisoBrittle use source_damage_anisoBrittle
use source_damage_anisoDuctile use source_damage_anisoDuctile
use results
implicit none implicit none
private private
integer, dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
damage_nonlocal_output damage_nonlocal_output
integer, dimension(:), allocatable, target, public :: & integer, dimension(:), allocatable, target, public :: &
@ -45,7 +44,7 @@ module damage_nonlocal
damage_nonlocal_getDiffusion33, & damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility, & damage_nonlocal_getMobility, &
damage_nonlocal_putNonLocalDamage, & damage_nonlocal_putNonLocalDamage, &
damage_nonlocal_postResults damage_nonlocal_Results
contains contains
@ -55,7 +54,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_init subroutine damage_nonlocal_init
integer :: maxNinstance,homog,instance,o,i integer :: maxNinstance,homog,instance,i
integer :: sizeState integer :: sizeState
integer :: NofMyHomog, h integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: & integer(kind(undefined_ID)) :: &
@ -69,7 +68,6 @@ subroutine damage_nonlocal_init
maxNinstance = count(damage_type == DAMAGE_nonlocal_ID) maxNinstance = count(damage_type == DAMAGE_nonlocal_ID)
if (maxNinstance == 0) return if (maxNinstance == 0) return
allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance)) allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance))
damage_nonlocal_output = '' damage_nonlocal_output = ''
allocate(damage_nonlocal_Noutput (maxNinstance), source=0) allocate(damage_nonlocal_Noutput (maxNinstance), source=0)
@ -92,7 +90,6 @@ subroutine damage_nonlocal_init
case ('damage') case ('damage')
damage_nonlocal_output(i,damage_typeInstance(h)) = outputs(i) damage_nonlocal_output(i,damage_typeInstance(h)) = outputs(i)
damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1 damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1
damage_nonlocal_sizePostResult(i,damage_typeInstance(h)) = 1
prm%outputID = [prm%outputID , damage_ID] prm%outputID = [prm%outputID , damage_ID]
end select end select
@ -107,7 +104,6 @@ subroutine damage_nonlocal_init
! allocate state arrays ! allocate state arrays
sizeState = 1 sizeState = 1
damageState(homog)%sizeState = sizeState damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog)) allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog)) allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog)) allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog))
@ -247,35 +243,29 @@ end subroutine damage_nonlocal_putNonLocalDamage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of damage results !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function damage_nonlocal_postResults(ip,el) subroutine damage_nonlocal_results(homog,group)
integer, intent(in) :: & integer, intent(in) :: homog
ip, & !< integration point character(len=*), intent(in) :: group
el !< element #if defined(PETSc) || defined(DAMASK_HDF5)
real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: & integer :: o, instance
damage_nonlocal_postResults
integer :: &
instance, homog, offset, o, c
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog) instance = damage_typeInstance(homog)
associate(prm => param(instance)) associate(prm => param(instance))
c = 0
outputsLoop: do o = 1,size(prm%outputID) outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o)) select case(prm%outputID(o))
case (damage_ID) case (damage_ID)
damage_nonlocal_postResults(c+1) = damage(homog)%p(offset) call results_writeDataset(group,damage(homog)%p,'phi',&
c = c + 1 'damage indicator','-')
end select end select
enddo outputsLoop enddo outputsLoop
end associate end associate
end function damage_nonlocal_postResults #endif
end subroutine damage_nonlocal_results
end module damage_nonlocal end module damage_nonlocal

78
src/grid/DAMASK_grid.f90
View File

@ -15,11 +15,7 @@ program DAMASK_spectral
use config use config
use debug use debug
use math use math
use mesh_grid
use CPFEM2 use CPFEM2
use FEsolving
use numerics
use homogenization
use material use material
use spectral_utilities use spectral_utilities
use grid_mech_spectral_basic use grid_mech_spectral_basic
@ -80,12 +76,6 @@ program DAMASK_spectral
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tLoadCase) :: newLoadCase type(tLoadCase) :: newLoadCase
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
integer(MPI_OFFSET_KIND) :: fileOffset
integer(MPI_OFFSET_KIND), dimension(:), allocatable :: outputSize
integer, parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742
integer, parameter :: maxRealOut = maxByteOut/pReal
integer(pLongInt), dimension(2) :: outputIndex
PetscErrorCode :: ierr
procedure(grid_mech_spectral_basic_init), pointer :: & procedure(grid_mech_spectral_basic_init), pointer :: &
mech_init mech_init
procedure(grid_mech_spectral_basic_forward), pointer :: & procedure(grid_mech_spectral_basic_forward), pointer :: &
@ -280,10 +270,8 @@ program DAMASK_spectral
enddo enddo
if (any(newLoadCase%stress%maskLogical .eqv. & if (any(newLoadCase%stress%maskLogical .eqv. &
newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only
if (any(newLoadCase%stress%maskLogical .and. & if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) &
transpose(newLoadCase%stress%maskLogical) .and. & .and. (math_I3<1))) errorID = 838 ! no rotation is allowed by stress BC
reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
errorID = 838 ! no rotation is allowed by stress BC
write(6,'(2x,a)') 'stress / GPa:' write(6,'(2x,a)') 'stress / GPa:'
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if(newLoadCase%stress%maskLogical(i,j)) then if(newLoadCase%stress%maskLogical(i,j)) then
@ -335,26 +323,10 @@ program DAMASK_spectral
! write header of output file ! write header of output file
if (worldrank == 0) then if (worldrank == 0) then
writeHeader: if (interface_restartInc < 1) then writeHeader: if (interface_restartInc < 1) then
open(newunit=fileUnit,file=trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED',status='REPLACE')
write(fileUnit) 'load:', trim(loadCaseFile) ! ... and write header
write(fileUnit) 'workingdir:', 'n/a'
write(fileUnit) 'geometry:', trim(geometryFile)
write(fileUnit) 'grid:', grid
write(fileUnit) 'size:', geomSize
write(fileUnit) 'materialpoint_sizeResults:', materialpoint_sizeResults
write(fileUnit) 'loadcases:', size(loadCases)
write(fileUnit) 'frequencies:', loadCases%outputfrequency ! one entry per LoadCase
write(fileUnit) 'times:', loadCases%time ! one entry per LoadCase
write(fileUnit) 'logscales:', loadCases%logscale
write(fileUnit) 'increments:', loadCases%incs ! one entry per LoadCase
write(fileUnit) 'startingIncrement:', interface_restartInc ! start with writing out the previous inc
write(fileUnit) 'eoh'
close(fileUnit) ! end of header
open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) & if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) &
write(6,'(/,a)') ' header of result and statistics file written out' write(6,'(/,a)') ' header of statistics file written out'
flush(6) flush(6)
else writeHeader else writeHeader
open(newunit=statUnit,file=trim(getSolverJobName())//& open(newunit=statUnit,file=trim(getSolverJobName())//&
@ -362,40 +334,11 @@ program DAMASK_spectral
endif writeHeader endif writeHeader
endif endif
!--------------------------------------------------------------------------------------------------
! prepare MPI parallel out (including opening of file)
allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND)
outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND)
call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_allreduce')
call MPI_file_open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.spectralOut', &
MPI_MODE_WRONLY + MPI_MODE_APPEND, &
MPI_INFO_NULL, &
fileUnit, &
ierr)
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_open')
call MPI_file_get_position(fileUnit,fileOffset,ierr) ! get offset from header
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_get_position')
fileOffset = fileOffset + sum(outputSize(1:worldrank)) ! offset of my process in file (header + processes before me)
call MPI_file_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr)
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_seek')
writeUndeformed: if (interface_restartInc < 1) then writeUndeformed: if (interface_restartInc < 1) then
write(6,'(1/,a)') ' ... writing initial configuration to file ........................' write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
call CPFEM_results(0,0.0_pReal) call CPFEM_results(0,0.0_pReal)
do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex = int([(i-1)*((maxRealOut)/materialpoint_sizeResults)+1, &
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
[(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)), &
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_write')
enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif writeUndeformed endif writeUndeformed
loadCaseLooping: do currentLoadCase = 1, size(loadCases) loadCaseLooping: do currentLoadCase = 1, size(loadCases)
time0 = time ! load case start time time0 = time ! load case start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc
@ -519,7 +462,6 @@ program DAMASK_spectral
write(6,'(/,a)') ' cutting back ' write(6,'(/,a)') ' cutting back '
else ! no more options to continue else ! no more options to continue
call IO_warning(850) call IO_warning(850)
call MPI_File_close(fileUnit,ierr)
close(statUnit) close(statUnit)
call quit(0) ! quit call quit(0) ! quit
endif endif
@ -537,19 +479,6 @@ program DAMASK_spectral
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................' write(6,'(1/,a)') ' ... writing results to file ......................................'
flush(6) flush(6)
call materialpoint_postResults()
call MPI_File_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr)
if (ierr /= 0) call IO_error(894, ext_msg='MPI_file_seek')
do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex=int([(i-1)*((maxRealOut)/materialpoint_sizeResults)+1, &
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),&
[(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)),&
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if(ierr /=0) call IO_error(894, ext_msg='MPI_file_write')
enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
call CPFEM_results(totalIncsCounter,time) call CPFEM_results(totalIncsCounter,time)
endif endif
if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then
@ -566,7 +495,6 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report summary of whole calculation ! report summary of whole calculation
write(6,'(/,a)') ' ###########################################################################' write(6,'(/,a)') ' ###########################################################################'
call MPI_file_close(fileUnit,ierr)
close(statUnit) close(statUnit)
call quit(0) ! no complains ;) call quit(0) ! no complains ;)

3
src/grid/grid_mech_FEM.f90
View File

@ -476,8 +476,7 @@ subroutine formResidual(da_local,x_local, &
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))

3
src/grid/grid_mech_spectral_basic.f90
View File

@ -440,8 +440,7 @@ subroutine formResidual(in, F, &
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))

3
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -509,8 +509,7 @@ subroutine formResidual(in, FandF_tau, &
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))

257
src/homogenization.f90
View File

@ -35,12 +35,6 @@ module homogenization
materialpoint_P !< first P--K stress of IP materialpoint_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: &
materialpoint_dPdF !< tangent of first P--K stress at IP materialpoint_dPdF !< tangent of first P--K stress at IP
real(pReal), dimension(:,:,:), allocatable, public :: &
materialpoint_results !< results array of material point
integer, public, protected :: &
materialpoint_sizeResults, &
thermal_maxSizePostResults, &
damage_maxSizePostResults
real(pReal), dimension(:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:), allocatable :: &
materialpoint_subF0, & !< def grad of IP at beginning of homogenization increment materialpoint_subF0, & !< def grad of IP at beginning of homogenization increment
@ -106,12 +100,12 @@ module homogenization
P,& !< partitioned stresses P,& !< partitioned stresses
F,& !< partitioned deformation gradients F,& !< partitioned deformation gradients
F0 !< partitioned initial deformation gradients F0 !< partitioned initial deformation gradients
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
real(pReal), dimension(3,3), intent(in) :: avgF !< average F real(pReal), dimension(3,3), intent(in) :: avgF !< average F
real(pReal), intent(in) :: dt !< time increment real(pReal), intent(in) :: dt !< time increment
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
end function mech_RGC_updateState end function mech_RGC_updateState
@ -125,7 +119,6 @@ module homogenization
public :: & public :: &
homogenization_init, & homogenization_init, &
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
materialpoint_postResults, &
homogenization_results homogenization_results
contains contains
@ -136,14 +129,6 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_init subroutine homogenization_init
integer, parameter :: FILEUNIT = 200
integer :: e,i,p
integer, dimension(:,:), pointer :: thisSize
integer, dimension(:) , pointer :: thisNoutput
character(len=64), dimension(:,:), pointer :: thisOutput
character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready
logical :: valid
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init
@ -156,80 +141,6 @@ subroutine homogenization_init
if (any(damage_type == DAMAGE_local_ID)) call damage_local_init if (any(damage_type == DAMAGE_local_ID)) call damage_local_init
if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init
!--------------------------------------------------------------------------------------------------
! write description file for homogenization output
mainProcess: if (worldrank == 0) then
call IO_write_jobFile(FILEUNIT,'outputHomogenization')
do p = 1,size(config_homogenization)
if (any(material_homogenizationAt == p)) then
write(FILEUNIT,'(/,a,/)') '['//trim(config_name_homogenization(p))//']'
write(FILEUNIT,'(a)') '(type) n/a'
write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p)
i = thermal_typeInstance(p) ! which instance of this thermal type
valid = .true. ! assume valid
select case(thermal_type(p)) ! split per thermal type
case (THERMAL_isothermal_ID)
outputName = THERMAL_isothermal_label
thisNoutput => null()
thisOutput => null()
thisSize => null()
case (THERMAL_adiabatic_ID)
outputName = THERMAL_adiabatic_label
thisNoutput => thermal_adiabatic_Noutput
thisOutput => thermal_adiabatic_output
thisSize => thermal_adiabatic_sizePostResult
case (THERMAL_conduction_ID)
outputName = THERMAL_conduction_label
thisNoutput => thermal_conduction_Noutput
thisOutput => thermal_conduction_output
thisSize => thermal_conduction_sizePostResult
case default
valid = .false.
end select
if (valid) then
write(FILEUNIT,'(a)') '(thermal)'//char(9)//trim(outputName)
if (thermal_type(p) /= THERMAL_isothermal_ID) then
do e = 1,thisNoutput(i)
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
endif
i = damage_typeInstance(p) ! which instance of this damage type
valid = .true. ! assume valid
select case(damage_type(p)) ! split per damage type
case (DAMAGE_none_ID)
outputName = DAMAGE_none_label
thisNoutput => null()
thisOutput => null()
thisSize => null()
case (DAMAGE_local_ID)
outputName = DAMAGE_local_label
thisNoutput => damage_local_Noutput
thisOutput => damage_local_output
thisSize => damage_local_sizePostResult
case (DAMAGE_nonlocal_ID)
outputName = DAMAGE_nonlocal_label
thisNoutput => damage_nonlocal_Noutput
thisOutput => damage_nonlocal_output
thisSize => damage_nonlocal_sizePostResult
case default
valid = .false.
end select
if (valid) then
write(FILEUNIT,'(a)') '(damage)'//char(9)//trim(outputName)
if (damage_type(p) /= DAMAGE_none_ID) then
do e = 1,thisNoutput(i)
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
endif
endif
enddo
close(FILEUNIT)
endif mainProcess
call config_deallocate('material.config/homogenization') call config_deallocate('material.config/homogenization')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -249,23 +160,7 @@ subroutine homogenization_init
allocate(materialpoint_converged(discretization_nIP,discretization_nElem), source=.true.) allocate(materialpoint_converged(discretization_nIP,discretization_nElem), source=.true.)
allocate(materialpoint_doneAndHappy(2,discretization_nIP,discretization_nElem), source=.true.) allocate(materialpoint_doneAndHappy(2,discretization_nIP,discretization_nElem), source=.true.)
!-------------------------------------------------------------------------------------------------- write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6)
! allocate and initialize global state and postresutls variables
thermal_maxSizePostResults = 0
damage_maxSizePostResults = 0
do p = 1,size(config_homogenization)
thermal_maxSizePostResults = max(thermal_maxSizePostResults, thermalState(p)%sizePostResults)
damage_maxSizePostResults = max(damage_maxSizePostResults, damageState (p)%sizePostResults)
enddo
materialpoint_sizeResults = 1 & ! grain count
+ 1 + thermal_maxSizePostResults &
+ damage_maxSizePostResults &
+ homogenization_maxNgrains * ( 1 & ! crystallite size
+ 1 + constitutive_source_maxSizePostResults)
allocate(materialpoint_results(materialpoint_sizeResults,discretization_nIP,discretization_nElem))
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
@ -581,52 +476,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
end subroutine materialpoint_stressAndItsTangent end subroutine materialpoint_stressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief parallelized calculation of result array at material points
!--------------------------------------------------------------------------------------------------
subroutine materialpoint_postResults
integer :: &
thePos, &
theSize, &
myNgrains, &
g, & !< grain number
i, & !< integration point number
e !< element number
!$OMP PARALLEL DO PRIVATE(myNgrains,thePos,theSize)
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
thePos = 0
theSize = thermalState (material_homogenizationAt(e))%sizePostResults &
+ damageState (material_homogenizationAt(e))%sizePostResults
materialpoint_results(thePos+1,i,e) = real(theSize,pReal) ! tell size of homogenization results
thePos = thePos + 1
if (theSize > 0) then ! any homogenization results to mention?
materialpoint_results(thePos+1:thePos+theSize,i,e) = postResults(i,e)
thePos = thePos + theSize
endif
materialpoint_results(thePos+1,i,e) = real(myNgrains,pReal) ! tell number of grains at materialpoint
thePos = thePos + 1
grainLooping :do g = 1,myNgrains
theSize = 1 + &
1 + &
sum(sourceState(material_phaseAt(g,e))%p(:)%sizePostResults)
materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(g,i,e) ! tell crystallite results
thePos = thePos + theSize
enddo grainLooping
enddo IpLooping
enddo elementLooping
!$OMP END PARALLEL DO
end subroutine materialpoint_postResults
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partition material point def grad onto constituents !> @brief partition material point def grad onto constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -738,81 +587,24 @@ subroutine averageStressAndItsTangent(ip,el)
end subroutine averageStressAndItsTangent end subroutine averageStressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion. call only,
!> if homogenization_sizePostResults(i,e) > 0 !!
!--------------------------------------------------------------------------------------------------
function postResults(ip,el)
integer, intent(in) :: &
ip, & !< integration point
el !< element number
real(pReal), dimension( thermalState (material_homogenizationAt(el))%sizePostResults &
+ damageState (material_homogenizationAt(el))%sizePostResults) :: &
postResults
integer :: &
startPos, endPos ,&
homog
postResults = 0.0_pReal
startPos = 1
endPos = thermalState(material_homogenizationAt(el))%sizePostResults
chosenThermal: select case (thermal_type(material_homogenizationAt(el)))
case (THERMAL_adiabatic_ID) chosenThermal
homog = material_homogenizationAt(el)
postResults(startPos:endPos) = &
thermal_adiabatic_postResults(homog,thermal_typeInstance(homog),thermalMapping(homog)%p(ip,el))
case (THERMAL_conduction_ID) chosenThermal
homog = material_homogenizationAt(el)
postResults(startPos:endPos) = &
thermal_conduction_postResults(homog,thermal_typeInstance(homog),thermalMapping(homog)%p(ip,el))
end select chosenThermal
startPos = endPos + 1
endPos = endPos + damageState(material_homogenizationAt(el))%sizePostResults
chosenDamage: select case (damage_type(material_homogenizationAt(el)))
case (DAMAGE_local_ID) chosenDamage
postResults(startPos:endPos) = damage_local_postResults(ip, el)
case (DAMAGE_nonlocal_ID) chosenDamage
postResults(startPos:endPos) = damage_nonlocal_postResults(ip, el)
end select chosenDamage
end function postResults
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes homogenization results to HDF5 output file !> @brief writes homogenization results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_results subroutine homogenization_results
use material, only: & use material, only: &
material_homogenization_type => homogenization_type material_homogenization_type => homogenization_type
integer :: p integer :: p
character(len=256) :: group character(len=pStringLen) :: group_base,group
!real(pReal), dimension(:,:,:), allocatable :: temp !real(pReal), dimension(:,:,:), allocatable :: temp
do p=1,size(config_name_homogenization) do p=1,size(config_name_homogenization)
group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p)) group_base = 'current/materialpoint/'//trim(config_name_homogenization(p))
call results_closeGroup(results_addGroup(group_base))
group = trim(group_base)//'/generic'
call results_closeGroup(results_addGroup(group)) call results_closeGroup(results_addGroup(group))
group = trim(group)//'/mech'
call results_closeGroup(results_addGroup(group))
select case(material_homogenization_type(p))
case(HOMOGENIZATION_rgc_ID)
call mech_RGC_results(homogenization_typeInstance(p),group)
end select
group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p))//'/generic'
call results_closeGroup(results_addGroup(group))
!temp = reshape(materialpoint_F,[3,3,discretization_nIP*discretization_nElem]) !temp = reshape(materialpoint_F,[3,3,discretization_nIP*discretization_nElem])
!call results_writeDataset(group,temp,'F',& !call results_writeDataset(group,temp,'F',&
! 'deformation gradient','1') ! 'deformation gradient','1')
@ -820,6 +612,31 @@ subroutine homogenization_results
!call results_writeDataset(group,temp,'P',& !call results_writeDataset(group,temp,'P',&
! '1st Piola-Kirchoff stress','Pa') ! '1st Piola-Kirchoff stress','Pa')
group = trim(group_base)//'/mech'
call results_closeGroup(results_addGroup(group))
select case(material_homogenization_type(p))
case(HOMOGENIZATION_rgc_ID)
call mech_RGC_results(homogenization_typeInstance(p),group)
end select
group = trim(group_base)//'/damage'
call results_closeGroup(results_addGroup(group))
select case(damage_type(p))
case(DAMAGE_LOCAL_ID)
call damage_local_results(p,group)
case(DAMAGE_NONLOCAL_ID)
call damage_nonlocal_results(p,group)
end select
group = trim(group_base)//'/thermal'
call results_closeGroup(results_addGroup(group))
select case(thermal_type(p))
case(THERMAL_ADIABATIC_ID)
call thermal_adiabatic_results(p,group)
case(THERMAL_CONDUCTION_ID)
call thermal_conduction_results(p,group)
end select
enddo enddo
end subroutine homogenization_results end subroutine homogenization_results

3
src/prec.f90
View File

@ -42,8 +42,7 @@ module prec
sizeState = 0, & !< size of state sizeState = 0, & !< size of state
sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
offsetDeltaState = 0, & !< index offset of delta state offsetDeltaState = 0, & !< index offset of delta state
sizeDeltaState = 0, & !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments sizeDeltaState = 0 !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments
sizePostResults = 0 !< size of output data
real(pReal), pointer, dimension(:), contiguous :: & real(pReal), pointer, dimension(:), contiguous :: &
atolState atolState
real(pReal), pointer, dimension(:,:), contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous real(pReal), pointer, dimension(:,:), contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous

52
src/source_damage_anisoBrittle.f90
View File

@ -13,6 +13,7 @@ module source_damage_anisoBrittle
use discretization use discretization
use config use config
use lattice use lattice
use results
implicit none implicit none
private private
@ -21,10 +22,7 @@ module source_damage_anisoBrittle
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
source_damage_anisoBrittle_instance !< instance of source mechanism source_damage_anisoBrittle_instance !< instance of source mechanism
integer, dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable :: &
source_damage_anisoBrittle_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_output !< name of each post result output source_damage_anisoBrittle_output !< name of each post result output
integer, dimension(:,:), allocatable :: & integer, dimension(:,:), allocatable :: &
@ -61,7 +59,7 @@ module source_damage_anisoBrittle
source_damage_anisoBrittle_init, & source_damage_anisoBrittle_init, &
source_damage_anisoBrittle_dotState, & source_damage_anisoBrittle_dotState, &
source_damage_anisobrittle_getRateAndItsTangent, & source_damage_anisobrittle_getRateAndItsTangent, &
source_damage_anisoBrittle_postResults source_damage_anisoBrittle_results
contains contains
@ -102,7 +100,6 @@ subroutine source_damage_anisoBrittle_init
enddo enddo
enddo enddo
allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0)
allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance)) allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoBrittle_output = '' source_damage_anisoBrittle_output = ''
@ -154,7 +151,6 @@ subroutine source_damage_anisoBrittle_init
select case(outputs(i)) select case(outputs(i))
case ('anisobrittle_drivingforce') case ('anisobrittle_drivingforce')
source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1
source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i) source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID] prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -171,7 +167,6 @@ subroutine source_damage_anisoBrittle_init
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
@ -262,39 +257,32 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
end subroutine source_damage_anisobrittle_getRateAndItsTangent end subroutine source_damage_anisoBrittle_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function source_damage_anisoBrittle_postResults(phase, constituent) subroutine source_damage_anisoBrittle_results(phase,group)
integer, intent(in) :: & integer, intent(in) :: phase
phase, & character(len=*), intent(in) :: group
constituent #if defined(PETSc) || defined(DAMASK_HDF5)
integer :: sourceOffset, o, instance
real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, & instance = source_damage_anisoBrittle_instance(phase)
source_damage_anisoBrittle_instance(phase)))) :: &
source_damage_anisoBrittle_postResults
integer :: &
instance, sourceOffset, o, c
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase) sourceOffset = source_damage_anisoBrittle_offset(phase)
c = 0 associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
outputsLoop: do o = 1,size(prm%outputID)
do o = 1,size(param(instance)%outputID) select case(prm%outputID(o))
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID) case (damage_drivingforce_ID)
source_damage_anisoBrittle_postResults(c+1) = & call results_writeDataset(group,stt,'tbd','driving force','tbd')
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
end select end select
enddo enddo outputsLoop
end function source_damage_anisoBrittle_postResults end associate
#endif
end subroutine source_damage_anisoBrittle_results
end module source_damage_anisoBrittle end module source_damage_anisoBrittle

48
src/source_damage_anisoDuctile.f90
View File

@ -12,6 +12,7 @@ module source_damage_anisoDuctile
use discretization use discretization
use material use material
use config use config
use results
implicit none implicit none
private private
@ -20,9 +21,6 @@ module source_damage_anisoDuctile
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism source_damage_anisoDuctile_instance !< instance of damage source mechanism
integer, dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_output !< name of each post result output source_damage_anisoDuctile_output !< name of each post result output
@ -54,7 +52,7 @@ module source_damage_anisoDuctile
source_damage_anisoDuctile_init, & source_damage_anisoDuctile_init, &
source_damage_anisoDuctile_dotState, & source_damage_anisoDuctile_dotState, &
source_damage_anisoDuctile_getRateAndItsTangent, & source_damage_anisoDuctile_getRateAndItsTangent, &
source_damage_anisoDuctile_postResults source_damage_anisoDuctile_results
contains contains
@ -96,7 +94,6 @@ subroutine source_damage_anisoDuctile_init
enddo enddo
enddo enddo
allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance)) allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoDuctile_output = '' source_damage_anisoDuctile_output = ''
@ -139,7 +136,6 @@ subroutine source_damage_anisoDuctile_init
select case(outputs(i)) select case(outputs(i))
case ('anisoductile_drivingforce') case ('anisoductile_drivingforce')
source_damage_anisoDuctile_sizePostResult(i,source_damage_anisoDuctile_instance(p)) = 1
source_damage_anisoDuctile_output(i,source_damage_anisoDuctile_instance(p)) = outputs(i) source_damage_anisoDuctile_output(i,source_damage_anisoDuctile_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID] prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -156,8 +152,7 @@ subroutine source_damage_anisoDuctile_init
sourceOffset = source_damage_anisoDuctile_offset(phase) sourceOffset = source_damage_anisoDuctile_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoDuctile_sizePostResult(:,instance)) sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
enddo enddo
@ -226,35 +221,28 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function source_damage_anisoDuctile_postResults(phase, constituent) subroutine source_damage_anisoDuctile_results(phase,group)
integer, intent(in) :: & integer, intent(in) :: phase
phase, & character(len=*), intent(in) :: group
constituent #if defined(PETSc) || defined(DAMASK_HDF5)
real(pReal), dimension(sum(source_damage_anisoDuctile_sizePostResult(:, & integer :: sourceOffset, o, instance
source_damage_anisoDuctile_instance(phase)))) :: &
source_damage_anisoDuctile_postResults
integer :: & instance = source_damage_anisoDuctile_instance(phase)
instance, sourceOffset, o, c
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase) sourceOffset = source_damage_anisoDuctile_offset(phase)
c = 0 associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
outputsLoop: do o = 1,size(prm%outputID)
do o = 1,size(param(instance)%outputID) select case(prm%outputID(o))
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID) case (damage_drivingforce_ID)
source_damage_anisoDuctile_postResults(c+1) = & call results_writeDataset(group,stt,'tbd','driving force','tbd')
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
end select end select
enddo enddo outputsLoop
end associate
#endif
end function source_damage_anisoDuctile_postResults end subroutine source_damage_anisoDuctile_results
end module source_damage_anisoDuctile end module source_damage_anisoDuctile

48
src/source_damage_isoBrittle.f90
View File

@ -12,15 +12,14 @@ module source_damage_isoBrittle
use discretization use discretization
use material use material
use config use config
use results
implicit none implicit none
private private
integer, dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
source_damage_isoBrittle_offset, & source_damage_isoBrittle_offset, &
source_damage_isoBrittle_instance source_damage_isoBrittle_instance
integer, dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable :: &
source_damage_isoBrittle_sizePostResult
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_isoBrittle_output source_damage_isoBrittle_output
enum, bind(c) enum, bind(c)
@ -46,7 +45,7 @@ module source_damage_isoBrittle
source_damage_isoBrittle_init, & source_damage_isoBrittle_init, &
source_damage_isoBrittle_deltaState, & source_damage_isoBrittle_deltaState, &
source_damage_isoBrittle_getRateAndItsTangent, & source_damage_isoBrittle_getRateAndItsTangent, &
source_damage_isoBrittle_postResults source_damage_isoBrittle_Results
contains contains
@ -86,7 +85,6 @@ subroutine source_damage_isoBrittle_init
enddo enddo
enddo enddo
allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),Ninstance)) allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),Ninstance))
source_damage_isoBrittle_output = '' source_damage_isoBrittle_output = ''
@ -122,7 +120,6 @@ subroutine source_damage_isoBrittle_init
select case(outputs(i)) select case(outputs(i))
case ('isobrittle_drivingforce') case ('isobrittle_drivingforce')
source_damage_isoBrittle_sizePostResult(i,source_damage_isoBrittle_instance(p)) = 1
source_damage_isoBrittle_output(i,source_damage_isoBrittle_instance(p)) = outputs(i) source_damage_isoBrittle_output(i,source_damage_isoBrittle_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID] prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -139,7 +136,6 @@ subroutine source_damage_isoBrittle_init
sourceOffset = source_damage_isoBrittle_offset(phase) sourceOffset = source_damage_isoBrittle_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,1) call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,1)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoBrittle_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
enddo enddo
@ -215,34 +211,30 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD
end subroutine source_damage_isoBrittle_getRateAndItsTangent end subroutine source_damage_isoBrittle_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function source_damage_isoBrittle_postResults(phase, constituent) subroutine source_damage_isoBrittle_results(phase,group)
integer, intent(in) :: & integer, intent(in) :: phase
phase, & character(len=*), intent(in) :: group
constituent #if defined(PETSc) || defined(DAMASK_HDF5)
real(pReal), dimension(sum(source_damage_isoBrittle_sizePostResult(:, & integer :: sourceOffset, o, instance
source_damage_isoBrittle_instance(phase)))) :: &
source_damage_isoBrittle_postResults
integer :: & instance = source_damage_isoBrittle_instance(phase)
instance, sourceOffset, o, c
instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase) sourceOffset = source_damage_isoBrittle_offset(phase)
c = 0 associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
outputsLoop: do o = 1,size(prm%outputID)
do o = 1,size(param(instance)%outputID) select case(prm%outputID(o))
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID) case (damage_drivingforce_ID)
source_damage_isoBrittle_postResults(c+1) = sourceState(phase)%p(sourceOffset)%state(1,constituent) call results_writeDataset(group,stt,'tbd','driving force','tbd')
c = c + 1
end select end select
enddo enddo outputsLoop
end function source_damage_isoBrittle_postResults end associate
#endif
end subroutine source_damage_isoBrittle_results
end module source_damage_isoBrittle end module source_damage_isoBrittle

51
src/source_damage_isoDuctile.f90
View File

@ -11,6 +11,7 @@ module source_damage_isoDuctile
use discretization use discretization
use material use material
use config use config
use results
implicit none implicit none
private private
@ -18,9 +19,6 @@ module source_damage_isoDuctile
source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_isoDuctile_instance !< instance of damage source mechanism source_damage_isoDuctile_instance !< instance of damage source mechanism
integer, dimension(:,:), allocatable, target, public :: &
source_damage_isoDuctile_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_isoDuctile_output !< name of each post result output source_damage_isoDuctile_output !< name of each post result output
@ -46,7 +44,7 @@ module source_damage_isoDuctile
source_damage_isoDuctile_init, & source_damage_isoDuctile_init, &
source_damage_isoDuctile_dotState, & source_damage_isoDuctile_dotState, &
source_damage_isoDuctile_getRateAndItsTangent, & source_damage_isoDuctile_getRateAndItsTangent, &
source_damage_isoDuctile_postResults source_damage_isoDuctile_Results
contains contains
@ -86,7 +84,6 @@ subroutine source_damage_isoDuctile_init
enddo enddo
enddo enddo
allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),Ninstance)) allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),Ninstance))
source_damage_isoDuctile_output = '' source_damage_isoDuctile_output = ''
@ -122,7 +119,6 @@ subroutine source_damage_isoDuctile_init
select case(outputs(i)) select case(outputs(i))
case ('isoductile_drivingforce') case ('isoductile_drivingforce')
source_damage_isoDuctile_sizePostResult(i,source_damage_isoDuctile_instance(p)) = 1
source_damage_isoDuctile_output(i,source_damage_isoDuctile_instance(p)) = outputs(i) source_damage_isoDuctile_output(i,source_damage_isoDuctile_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID] prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -138,10 +134,8 @@ subroutine source_damage_isoDuctile_init
sourceOffset = source_damage_isoDuctile_offset(phase) sourceOffset = source_damage_isoDuctile_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoDuctile_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
enddo enddo
end subroutine source_damage_isoDuctile_init end subroutine source_damage_isoDuctile_init
@ -197,34 +191,31 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD
end subroutine source_damage_isoDuctile_getRateAndItsTangent end subroutine source_damage_isoDuctile_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function source_damage_isoDuctile_postResults(phase, constituent) subroutine source_damage_isoDuctile_results(phase,group)
integer, intent(in) :: & integer, intent(in) :: phase
phase, & character(len=*), intent(in) :: group
constituent #if defined(PETSc) || defined(DAMASK_HDF5)
real(pReal), dimension(sum(source_damage_isoDuctile_sizePostResult(:, & integer :: sourceOffset, o, instance
source_damage_isoDuctile_instance(phase)))) :: &
source_damage_isoDuctile_postResults
integer :: & instance = source_damage_isoDuctile_instance(phase)
instance, sourceOffset, o, c sourceOffset = source_damage_isoDuctile_offset(phase)
instance = source_damage_isoDuctile_instance(phase) associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
sourceOffset = source_damage_isoDuctile_offset(phase) outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_drivingforce_ID)
call results_writeDataset(group,stt,'tbd','driving force','tbd')
end select
enddo outputsLoop
end associate
#endif
c = 0 end subroutine source_damage_isoDuctile_results
do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID)
source_damage_isoDuctile_postResults(c+1) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
end select
enddo
end function source_damage_isoDuctile_postResults
end module source_damage_isoDuctile end module source_damage_isoDuctile

43
src/thermal_adiabatic.f90
View File

@ -7,6 +7,7 @@ module thermal_adiabatic
use config use config
use numerics use numerics
use material use material
use results
use source_thermal_dissipation use source_thermal_dissipation
use source_thermal_externalheat use source_thermal_externalheat
use crystallite use crystallite
@ -15,8 +16,6 @@ module thermal_adiabatic
implicit none implicit none
private private
integer, dimension(:,:), allocatable, target, public :: &
thermal_adiabatic_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
thermal_adiabatic_output !< name of each post result output thermal_adiabatic_output !< name of each post result output
@ -37,7 +36,7 @@ module thermal_adiabatic
thermal_adiabatic_getSourceAndItsTangent, & thermal_adiabatic_getSourceAndItsTangent, &
thermal_adiabatic_getSpecificHeat, & thermal_adiabatic_getSpecificHeat, &
thermal_adiabatic_getMassDensity, & thermal_adiabatic_getMassDensity, &
thermal_adiabatic_postResults thermal_adiabatic_results
contains contains
@ -57,7 +56,6 @@ subroutine thermal_adiabatic_init
maxNinstance = count(thermal_type == THERMAL_adiabatic_ID) maxNinstance = count(thermal_type == THERMAL_adiabatic_ID)
if (maxNinstance == 0) return if (maxNinstance == 0) return
allocate(thermal_adiabatic_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(thermal_adiabatic_output (maxval(homogenization_Noutput),maxNinstance)) allocate(thermal_adiabatic_output (maxval(homogenization_Noutput),maxNinstance))
thermal_adiabatic_output = '' thermal_adiabatic_output = ''
allocate(thermal_adiabatic_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) allocate(thermal_adiabatic_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
@ -75,14 +73,12 @@ subroutine thermal_adiabatic_init
thermal_adiabatic_Noutput(instance) = thermal_adiabatic_Noutput(instance) + 1 thermal_adiabatic_Noutput(instance) = thermal_adiabatic_Noutput(instance) + 1
thermal_adiabatic_outputID(thermal_adiabatic_Noutput(instance),instance) = temperature_ID thermal_adiabatic_outputID(thermal_adiabatic_Noutput(instance),instance) = temperature_ID
thermal_adiabatic_output(thermal_adiabatic_Noutput(instance),instance) = outputs(i) thermal_adiabatic_output(thermal_adiabatic_Noutput(instance),instance) = outputs(i)
thermal_adiabatic_sizePostResult(thermal_adiabatic_Noutput(instance),instance) = 1
end select end select
enddo enddo
! allocate state arrays ! allocate state arrays
sizeState = 1 sizeState = 1
thermalState(section)%sizeState = sizeState thermalState(section)%sizeState = sizeState
thermalState(section)%sizePostResults = sum(thermal_adiabatic_sizePostResult(:,instance))
allocate(thermalState(section)%state0 (sizeState,NofMyHomog), source=thermal_initialT(section)) allocate(thermalState(section)%state0 (sizeState,NofMyHomog), source=thermal_initialT(section))
allocate(thermalState(section)%subState0(sizeState,NofMyHomog), source=thermal_initialT(section)) allocate(thermalState(section)%subState0(sizeState,NofMyHomog), source=thermal_initialT(section))
allocate(thermalState(section)%state (sizeState,NofMyHomog), source=thermal_initialT(section)) allocate(thermalState(section)%state (sizeState,NofMyHomog), source=thermal_initialT(section))
@ -252,32 +248,27 @@ end function thermal_adiabatic_getMassDensity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of thermal results !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function thermal_adiabatic_postResults(homog,instance,of) result(postResults) subroutine thermal_adiabatic_results(homog,group)
integer, intent(in) :: & integer, intent(in) :: homog
homog, & character(len=*), intent(in) :: group
instance, & #if defined(PETSc) || defined(DAMASK_HDF5)
of integer :: o, instance
real(pReal), dimension(sum(thermal_adiabatic_sizePostResult(:,instance))) :: & instance = thermal_typeInstance(homog)
postResults
integer :: & outputsLoop: do o = 1,thermal_adiabatic_Noutput(instance)
o, c
c = 0
do o = 1,thermal_adiabatic_Noutput(instance)
select case(thermal_adiabatic_outputID(o,instance)) select case(thermal_adiabatic_outputID(o,instance))
case (temperature_ID) case (temperature_ID)
postResults(c+1) = temperature(homog)%p(of) call results_writeDataset(group,temperature(homog)%p,'T',&
c = c + 1 'temperature','K')
end select end select
enddo enddo outputsLoop
#endif
end function thermal_adiabatic_postResults end subroutine thermal_adiabatic_results
end module thermal_adiabatic end module thermal_adiabatic

42
src/thermal_conduction.f90
View File

@ -7,6 +7,7 @@ module thermal_conduction
use material use material
use config use config
use lattice use lattice
use results
use crystallite use crystallite
use source_thermal_dissipation use source_thermal_dissipation
use source_thermal_externalheat use source_thermal_externalheat
@ -14,8 +15,6 @@ module thermal_conduction
implicit none implicit none
private private
integer, dimension(:,:), allocatable, target, public :: &
thermal_conduction_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
thermal_conduction_output !< name of each post result output thermal_conduction_output !< name of each post result output
@ -37,7 +36,7 @@ module thermal_conduction
thermal_conduction_getSpecificHeat, & thermal_conduction_getSpecificHeat, &
thermal_conduction_getMassDensity, & thermal_conduction_getMassDensity, &
thermal_conduction_putTemperatureAndItsRate, & thermal_conduction_putTemperatureAndItsRate, &
thermal_conduction_postResults thermal_conduction_results
contains contains
@ -60,7 +59,6 @@ subroutine thermal_conduction_init
maxNinstance = count(thermal_type == THERMAL_conduction_ID) maxNinstance = count(thermal_type == THERMAL_conduction_ID)
if (maxNinstance == 0) return if (maxNinstance == 0) return
allocate(thermal_conduction_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(thermal_conduction_output (maxval(homogenization_Noutput),maxNinstance)) allocate(thermal_conduction_output (maxval(homogenization_Noutput),maxNinstance))
thermal_conduction_output = '' thermal_conduction_output = ''
allocate(thermal_conduction_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) allocate(thermal_conduction_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
@ -78,7 +76,6 @@ subroutine thermal_conduction_init
thermal_conduction_Noutput(instance) = thermal_conduction_Noutput(instance) + 1 thermal_conduction_Noutput(instance) = thermal_conduction_Noutput(instance) + 1
thermal_conduction_outputID(thermal_conduction_Noutput(instance),instance) = temperature_ID thermal_conduction_outputID(thermal_conduction_Noutput(instance),instance) = temperature_ID
thermal_conduction_output(thermal_conduction_Noutput(instance),instance) = outputs(i) thermal_conduction_output(thermal_conduction_Noutput(instance),instance) = outputs(i)
thermal_conduction_sizePostResult(thermal_conduction_Noutput(instance),instance) = 1
end select end select
enddo enddo
@ -86,7 +83,6 @@ subroutine thermal_conduction_init
! allocate state arrays ! allocate state arrays
sizeState = 0 sizeState = 0
thermalState(section)%sizeState = sizeState thermalState(section)%sizeState = sizeState
thermalState(section)%sizePostResults = sum(thermal_conduction_sizePostResult(:,instance))
allocate(thermalState(section)%state0 (sizeState,NofMyHomog)) allocate(thermalState(section)%state0 (sizeState,NofMyHomog))
allocate(thermalState(section)%subState0(sizeState,NofMyHomog)) allocate(thermalState(section)%subState0(sizeState,NofMyHomog))
allocate(thermalState(section)%state (sizeState,NofMyHomog)) allocate(thermalState(section)%state (sizeState,NofMyHomog))
@ -265,31 +261,27 @@ end subroutine thermal_conduction_putTemperatureAndItsRate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of thermal results !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function thermal_conduction_postResults(homog,instance,of) result(postResults) subroutine thermal_conduction_results(homog,group)
integer, intent(in) :: & integer, intent(in) :: homog
homog, & character(len=*), intent(in) :: group
instance, & #if defined(PETSc) || defined(DAMASK_HDF5)
of integer :: o, instance
real(pReal), dimension(sum(thermal_conduction_sizePostResult(:,instance))) :: & instance = thermal_typeInstance(homog)
postResults
integer :: & outputsLoop: do o = 1,thermal_conduction_Noutput(instance)
o, c
c = 0
do o = 1,thermal_conduction_Noutput(instance)
select case(thermal_conduction_outputID(o,instance)) select case(thermal_conduction_outputID(o,instance))
case (temperature_ID) case (temperature_ID)
postResults(c+1) = temperature(homog)%p(of) call results_writeDataset(group,temperature(homog)%p,'T',&
c = c + 1 'temperature','K')
end select end select
enddo enddo outputsLoop
#endif
end function thermal_conduction_postResults end subroutine thermal_conduction_results
end module thermal_conduction end module thermal_conduction