Merge branch 'development' into even-more-HDF5-postprocessing

This commit is contained in:
Martin Diehl 2019-07-06 20:35:09 -07:00
commit 18f00cb8df
104 changed files with 4221 additions and 4186 deletions

View File

@ -110,7 +110,7 @@ for executable in icc icpc ifort ;do
done
firstLevel "MPI Wrappers"
for executable in mpicc mpiCC mpic++ mpicpc mpicxx mpifort mpif90 mpif77; do
for executable in mpicc mpiCC mpiicc mpic++ mpicpc mpicxx mpifort mpif90 mpif77; do
getDetails $executable '-show'
done

@ -1 +1 @@
Subproject commit 1d3cf8180a20bcba6958ce82eb97befec077d7d2
Subproject commit 18a976753be06aca6e15f580998e713daa08bb41

View File

@ -1 +1 @@
v2.0.3-559-g589f5343
v2.0.3-614-g4d6a047b

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,12 +1,19 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,re,sys,collections
import math,scipy,scipy.linalg # noqa
import numpy as np
import os
import sys
from optparse import OptionParser
import re
import collections
import math # noqa
import scipy # noqa
import scipy.linalg # noqa
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,13 +1,15 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os
import math
from optparse import OptionParser
import numpy as np
import scipy.ndimage
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -27,7 +30,7 @@ def curlFFT(geomdim,field):
curl_fourier = np.empty(field_fourier.shape,'c16')
# differentiation in Fourier space
TWOPIIMG = 2.0j*math.pi
TWOPIIMG = 2.0j*np.pi
einsums = {
3:'slm,ijkl,ijkm->ijks', # vector, 3 -> 3
9:'slm,ijkl,ijknm->ijksn', # tensor, 3x3 -> 3x3

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -30,6 +33,7 @@ def derivative(coordinates,what):
return result
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,10 +1,12 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import os
import sys
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -16,6 +18,7 @@ def determinant(m):
-m[1]*m[3]*m[8] \
-m[0]*m[5]*m[7]
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,8 +1,9 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import os
import sys
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
@ -19,6 +20,7 @@ def deviator(m,spherical = False):
]
return dev,sph if spherical else dev
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,12 +1,15 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import os
import sys
from optparse import OptionParser
import numpy as np
import scipy.ndimage
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -57,7 +60,7 @@ def displacementAvgFFT(F,grid,size,nodal=False,transformed=False):
#--------------------------------------------------------------------------------------------------
def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
"""Calculate cell center (or nodal) displacement for deformation gradient field specified in each grid cell"""
integrator = 0.5j * size / math.pi
integrator = 0.5j * size / np.pi
kk, kj, ki = np.meshgrid(np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2])),
np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1])),

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -27,7 +30,7 @@ def divFFT(geomdim,field):
div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16')
# differentiation in Fourier space
TWOPIIMG = 2.0j*math.pi
TWOPIIMG = 2.0j*np.pi
einsums = {
3:'ijkl,ijkl->ijk', # vector, 3 -> 1
9:'ijkm,ijklm->ijkl', # tensor, 3x3 -> 3

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -26,6 +29,7 @@ def E_hkl(stiffness,vec): # stiffness = (c11,c12,c44)
return 1.0/invE
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,12 +1,16 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,itertools
import os
import sys
from optparse import OptionParser
import itertools
import numpy as np
from scipy import ndimage
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -31,6 +35,7 @@ def periodic_3Dpad(array, rimdim=(1,1,1)):
padded[p[0],p[1],p[2]] = array[spot[0],spot[1],spot[2]]
return padded
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,12 +1,15 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
from scipy import ndimage
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -27,7 +30,7 @@ def gradFFT(geomdim,field):
grad_fourier = np.empty(field_fourier.shape+(3,),'c16')
# differentiation in Fourier space
TWOPIIMG = 2.0j*math.pi
TWOPIIMG = 2.0j*np.pi
einsums = {
1:'ijkl,ijkm->ijkm', # scalar, 1 -> 3
3:'ijkl,ijkm->ijklm', # vector, 3 -> 3x3

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,10 +1,11 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,12 +1,15 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
import os
import sys
from optparse import OptionParser
from collections import OrderedDict
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -14,7 +17,7 @@ def Mises(what,tensor):
dev = tensor - np.trace(tensor)/3.0*np.eye(3)
symdev = 0.5*(dev+dev.T)
return math.sqrt(np.sum(symdev*symdev.T)*
return np.sqrt(np.sum(symdev*symdev.T)*
{
'stress': 3.0/2.0,
'strain': 2.0/3.0,

View File

@ -1,10 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -16,12 +20,13 @@ def norm(which,object):
if which == 'Abs': # p = 1
return sum(map(abs, object))
elif which == 'Frobenius': # p = 2
return math.sqrt(sum([x*x for x in object]))
return np.sqrt(sum([x*x for x in object]))
elif which == 'Max': # p = inf
return max(map(abs, object))
else:
return -1
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
@ -136,7 +139,7 @@ parser.set_defaults(force = (0.0,0.0,1.0),
quaternion='orientation',
normal = None,
lattice = latticeChoices[0],
CoverA = math.sqrt(8./3.),
CoverA = np.sqrt(8./3.),
)
(options, filenames) = parser.parse_args()

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,19 +1,21 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import os
import sys
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
Append data of ASCIItable(s).
Append data of ASCIItable(s) column-wise.
""", version = scriptID)

View File

@ -1,15 +1,19 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import os
import sys
from optparse import OptionParser
import numpy as np
import scipy.ndimage
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,12 +1,17 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,re,sys,fnmatch
import math # noqa
import numpy as np
import os
import sys
from optparse import OptionParser
import re
import fnmatch
import math # noqa
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,12 +1,15 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import math # noqa
import numpy as np
import os
import sys
from optparse import OptionParser, OptionGroup
import math # noqa
import numpy as np
import damask
def periodicAverage(coords, limits):
"""Centroid in periodic domain, see https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions"""
theta = 2.0*np.pi * (coords - limits[0])/(limits[1] - limits[0])

79
processing/post/growTable.py Executable file
View File

@ -0,0 +1,79 @@
#!/usr/bin/env python3
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
Append data of ASCIItable(s) row-wise.
""", version = scriptID)
parser.add_option('-a', '--add','--table',
dest = 'table',
action = 'extend', metavar = '<string LIST>',
help = 'tables to add')
(options,filenames) = parser.parse_args()
if options.table is None:
parser.error('no table specified.')
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try: table = damask.ASCIItable(name = name,
buffered = False)
except: continue
damask.util.report(scriptName,name)
tables = []
for addTable in options.table:
try: tables.append(damask.ASCIItable(name = addTable,
buffered = False,
readonly = True)
)
except: continue
# ------------------------------------------ read headers ------------------------------------------
table.head_read()
for addTable in tables: addTable.head_read()
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.head_write()
# ------------------------------------------ process data ------------------------------------------
table.data_readArray()
data = table.data
for addTable in tables:
addTable.data_readArray(table.labels(raw = True))
data = np.vstack((data,addTable.data))
table.data = data
table.data_writeArray()
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables
for addTable in tables:
addTable.close()

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,15 +1,19 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
from PIL import Image
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,15 +1,19 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
from PIL import Image, ImageDraw
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,15 +1,19 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
from PIL import Image
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,10 +1,13 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys
import damask
import os
import sys
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,13 +1,16 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import sys,os
import damask
import os
import sys
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,9 +1,15 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math,re,time,struct
import damask
import os
import sys
from optparse import OptionParser, OptionGroup
import math
import re
import time
import struct
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,13 +1,17 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,re
import damask
import os
import sys
from optparse import OptionParser
import re
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,14 +1,18 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,13 +1,15 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,12 +1,17 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,string,math,sys
import numpy as np
import os
import sys
from optparse import OptionParser
import string
import math
import numpy as np
import vtk
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])

View File

@ -1,15 +1,19 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,vtk
import damask
from vtk.util import numpy_support
from collections import defaultdict
import os
from optparse import OptionParser
from collections import defaultdict
import vtk
from vtk.util import numpy_support
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,15 +1,19 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,vtk
import damask
from collections import defaultdict
import os
from optparse import OptionParser
from collections import defaultdict
import vtk
from vtk.util import numpy_support
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,15 +1,20 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,vtk
import damask
from vtk.util import numpy_support
from collections import defaultdict
import os
from optparse import OptionParser
from collections import defaultdict
import vtk
from vtk.util import numpy_support
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,14 +1,19 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,vtk
import numpy as np
import damask
import os
import sys
from optparse import OptionParser
import vtk
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -1,14 +1,19 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,vtk
import numpy as np
import damask
import os
import sys
from optparse import OptionParser
import vtk
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------

View File

@ -125,11 +125,6 @@ subroutine CPFEM_init
! flush(6)
! endif
! call IO_read_intFile(777,'recordedPhase'//trim(rankStr),modelName,size(material_phase))
! read (777,rec=1) material_phase
! close (777)
! call IO_read_realFile(777,'convergedF'//trim(rankStr),modelName,size(crystallite_F0))
! read (777,rec=1) crystallite_F0
! close (777)
@ -262,7 +257,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
if (debug_e <= discretization_nElem .and. debug_i <=discretization_nIP) then
write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)),/)') &
'<< CPFEM >> aged state of elFE ip grain',debug_e, debug_i, 1, &
plasticState(phaseAt(1,debug_i,debug_e))%state(:,phasememberAt(1,debug_i,debug_e))
plasticState(material_phaseAt(1,debug_e))%state(:,material_phasememberAt(1,debug_i,debug_e))
endif
endif
@ -280,10 +275,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
! write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files'
!
! call IO_write_jobRealFile(777,'recordedPhase'//trim(rankStr),size(material_phase))
! write (777,rec=1) material_phase
! close (777)
! call IO_write_jobRealFile(777,'convergedF'//trim(rankStr),size(crystallite_F0))
! write (777,rec=1) crystallite_F0
! close (777)

View File

@ -42,7 +42,7 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll
call DAMASK_interface_init ! Spectral and FEM interface to commandline
call DAMASK_interface_init ! Spectral and FEM interface to commandline
call prec_init
call IO_init
#ifdef FEM
@ -52,7 +52,6 @@ subroutine CPFEM_initAll
call debug_init
call config_init
call math_init
call FE_init
call mesh_init
call lattice_init
call HDF5_utilities_init
@ -78,8 +77,8 @@ subroutine CPFEM_init
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
flush(6)
! *** restore the last converged values of each essential variable from the binary file
if (restartRead) then
! *** restore the last converged values of each essential variable
if (interface_restartInc > 0) then
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0) then
write(6,'(a)') '<< CPFEM >> restored state variables of last converged step from hdf5 file'
flush(6)
@ -89,31 +88,28 @@ subroutine CPFEM_init
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
call HDF5_read(fileHandle,material_phase, 'recordedPhase')
call HDF5_read(fileHandle,crystallite_F0, 'convergedF')
call HDF5_read(fileHandle,crystallite_Fp0, 'convergedFp')
call HDF5_read(fileHandle,crystallite_Fi0, 'convergedFi')
call HDF5_read(fileHandle,crystallite_Lp0, 'convergedLp')
call HDF5_read(fileHandle,crystallite_Li0, 'convergedLi')
call HDF5_read(fileHandle,crystallite_S0, 'convergedS')
call HDF5_read(fileHandle,crystallite_F0, 'convergedF')
call HDF5_read(fileHandle,crystallite_Fp0,'convergedFp')
call HDF5_read(fileHandle,crystallite_Fi0,'convergedFi')
call HDF5_read(fileHandle,crystallite_Lp0,'convergedLp')
call HDF5_read(fileHandle,crystallite_Li0,'convergedLi')
call HDF5_read(fileHandle,crystallite_S0, 'convergedS')
groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases')
do ph = 1,size(phase_plasticity)
write(PlasticItem,*) ph,'_'
call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst')
write(PlasticItem,*) ph,'_'
call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst')
enddo
call HDF5_closeGroup(groupPlasticID)
groupHomogID = HDF5_openGroup(fileHandle,'HomogStates')
do homog = 1, material_Nhomogenization
write(HomogItem,*) homog,'_'
call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'convergedStateHomog')
write(HomogItem,*) homog,'_'
call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'convergedStateHomog')
enddo
call HDF5_closeGroup(groupHomogID)
call HDF5_closeFile(fileHandle)
restartRead = .false.
endif
end subroutine CPFEM_init
@ -136,7 +132,7 @@ subroutine CPFEM_age
crystallite_Lp0 = crystallite_Lp
crystallite_Fi0 = crystallite_Fi
crystallite_Li0 = crystallite_Li
crystallite_S0 = crystallite_S
crystallite_S0 = crystallite_S
do i = 1, size(plasticState)
plasticState(i)%state0 = plasticState(i)%state
@ -158,13 +154,12 @@ subroutine CPFEM_age
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a')
call HDF5_write(fileHandle,material_phase, 'recordedPhase')
call HDF5_write(fileHandle,crystallite_F0, 'convergedF')
call HDF5_write(fileHandle,crystallite_Fp0, 'convergedFp')
call HDF5_write(fileHandle,crystallite_Fi0, 'convergedFi')
call HDF5_write(fileHandle,crystallite_Lp0, 'convergedLp')
call HDF5_write(fileHandle,crystallite_Li0, 'convergedLi')
call HDF5_write(fileHandle,crystallite_S0, 'convergedS')
call HDF5_write(fileHandle,crystallite_F0, 'convergedF')
call HDF5_write(fileHandle,crystallite_Fp0, 'convergedFp')
call HDF5_write(fileHandle,crystallite_Fi0, 'convergedFi')
call HDF5_write(fileHandle,crystallite_Lp0, 'convergedLp')
call HDF5_write(fileHandle,crystallite_Li0, 'convergedLi')
call HDF5_write(fileHandle,crystallite_S0, 'convergedS')
groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases')
do ph = 1,size(phase_plasticity)
@ -195,8 +190,8 @@ end subroutine CPFEM_age
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_results(inc,time)
integer, intent(in) :: inc
real(pReal), intent(in) :: time
integer, intent(in) :: inc
real(pReal), intent(in) :: time
call results_openJobFile
call results_addIncrement(inc,time)

View File

@ -353,8 +353,7 @@ subroutine flux(f,ts,n,time)
!--------------------------------------------------------------------------------------------------
!> @brief sets user defined output variables for Marc
!> @details select a variable contour plotting (user subroutine).
!> @brief trigger writing of results
!--------------------------------------------------------------------------------------------------
subroutine uedinc(inc,incsub)
use prec

View File

@ -1,8 +1,7 @@
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief triggering reading in of restart information when doing a restart
!> @todo Descriptions for public variables needed
!> @brief holds some global variables and gets extra information for commercial FEM
!--------------------------------------------------------------------------------------------------
module FEsolving
use prec
@ -12,32 +11,34 @@ module FEsolving
implicit none
private
integer, public :: &
restartInc = 1 !< needs description
logical, public :: &
symmetricSolver = .false., & !< use a symmetric FEM solver
restartWrite = .false., & !< write current state to enable restart
#if defined(Marc4DAMASK) || defined(Abaqus)
restartRead = .false., & !< restart information to continue calculation from saved state
#endif
restartWrite = .false., & !< write current state to enable restart
terminallyIll = .false. !< at least one material point is terminally ill
integer, dimension(:,:), allocatable, public :: &
FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP
integer, dimension(2), public :: &
integer, dimension(2), public :: &
FEsolving_execElem !< for ping-pong scheme always whole range, otherwise one specific element
#if defined(Marc4DAMASK) || defined(Abaqus)
logical, public, protected :: &
symmetricSolver = .false. !< use a symmetric FEM solver (only Abaqus)
character(len=1024), public :: &
modelName !< needs description
logical, dimension(:,:), allocatable, public :: &
calcMode !< do calculation or simply collect when using ping pong scheme
public :: FE_init
#endif
contains
#if defined(Marc4DAMASK) || defined(Abaqus)
!--------------------------------------------------------------------------------------------------
!> @brief determine whether a symmetric solver is used and whether restart is requested
!> @details restart information is found in input file in case of FEM solvers, in case of spectal
@ -45,27 +46,15 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine FE_init
#if defined(Marc4DAMASK) || defined(Abaqus)
integer, parameter :: &
FILEUNIT = 222
integer :: j
character(len=65536) :: tag, line
integer, allocatable, dimension(:) :: chunkPos
#endif
write(6,'(/,a)') ' <<<+- FEsolving init -+>>>'
modelName = getSolverJobName()
#if defined(Grid) || defined(FEM)
restartInc = interface_RestartInc
if(restartInc < 0) then
call IO_warning(warning_ID=34)
restartInc = 0
endif
restartRead = restartInc > 0 ! only read in if "true" restart requested
#else
call IO_open_inputFile(FILEUNIT,modelName)
rewind(FILEUNIT)
do
@ -125,7 +114,6 @@ subroutine FE_init
200 close(FILEUNIT)
endif
#endif
if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0) then
write(6,'(a21,l1)') ' restart writing: ', restartWrite
write(6,'(a21,l1)') ' restart reading: ', restartRead
@ -133,5 +121,6 @@ subroutine FE_init
endif
end subroutine FE_init
#endif
end module FEsolving

View File

@ -22,7 +22,7 @@ module HDF5_utilities
#if defined(PETSc) || defined(DAMASK_HDF5)
!--------------------------------------------------------------------------------------------------
!> @brief reads integer or float data of defined shape from file ! ToDo: order of arguments wrong
!> @brief reads integer or float data of defined shape from file ! ToDo: order of arguments wrong
!> @details for parallel IO, all dimension except for the last need to match
!--------------------------------------------------------------------------------------------------
interface HDF5_read
@ -45,7 +45,7 @@ module HDF5_utilities
end interface HDF5_read
!--------------------------------------------------------------------------------------------------
!> @brief writes integer or real data of defined shape to file ! ToDo: order of arguments wrong
!> @brief writes integer or real data of defined shape to file ! ToDo: order of arguments wrong
!> @details for parallel IO, all dimension except for the last need to match
!--------------------------------------------------------------------------------------------------
interface HDF5_write
@ -1759,66 +1759,66 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
myStart, globalShape, &
loc_id,localShape,datasetName,parallel)
integer(HID_T), intent(in) :: loc_id !< file or group handle
character(len=*), intent(in) :: datasetName !< name of the dataset in the file
logical, intent(in) :: parallel
integer(HSIZE_T), intent(in), dimension(:) :: &
localShape
integer(HSIZE_T), intent(out), dimension(size(localShape,1)):: &
myStart, &
globalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
integer(HID_T), intent(in) :: loc_id !< file or group handle
character(len=*), intent(in) :: datasetName !< name of the dataset in the file
logical, intent(in) :: parallel
integer(HSIZE_T), intent(in), dimension(:) :: &
localShape
integer(HSIZE_T), intent(out), dimension(size(localShape,1)):: &
myStart, &
globalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
integer, dimension(worldsize) :: &
readSize !< contribution of all processes
integer :: ierr
integer :: hdferr
integer, dimension(worldsize) :: &
readSize !< contribution of all processes
integer :: ierr
integer :: hdferr
!-------------------------------------------------------------------------------------------------
! creating a property list for transfer properties (is collective for MPI)
call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pcreate_f')
call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pcreate_f')
!--------------------------------------------------------------------------------------------------
readSize = 0
readSize(worldrank+1) = int(localShape(ubound(localShape,1)))
readSize = 0
readSize(worldrank+1) = int(localShape(ubound(localShape,1)))
#ifdef PETSc
if (parallel) then
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pset_dxpl_mpio_f')
call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if (ierr /= 0) call IO_error(894,ext_msg='initialize_read: MPI_allreduce')
endif
if (parallel) then
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pset_dxpl_mpio_f')
call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if (ierr /= 0) call IO_error(894,ext_msg='initialize_read: MPI_allreduce')
endif
#endif
myStart = int(0,HSIZE_T)
myStart(ubound(myStart)) = int(sum(readSize(1:worldrank)),HSIZE_T)
globalShape = [localShape(1:ubound(localShape,1)-1),int(sum(readSize),HSIZE_T)]
myStart = int(0,HSIZE_T)
myStart(ubound(myStart)) = int(sum(readSize(1:worldrank)),HSIZE_T)
globalShape = [localShape(1:ubound(localShape,1)-1),int(sum(readSize),HSIZE_T)]
!--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape)
call h5screate_simple_f(size(localShape), localShape, memspace_id, hdferr, localShape)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5screate_simple_f/memspace_id')
call h5screate_simple_f(size(localShape), localShape, memspace_id, hdferr, localShape)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5screate_simple_f/memspace_id')
!--------------------------------------------------------------------------------------------------
! creating a property list for IO and set it to collective
call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pcreate_f')
call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pcreate_f')
#ifdef PETSc
call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pset_all_coll_metadata_ops_f')
call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pset_all_coll_metadata_ops_f')
#endif
!--------------------------------------------------------------------------------------------------
! open the dataset in the file and get the space ID
call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5dopen_f')
call h5dget_space_f(dset_id, filespace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5dget_space_f')
call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5dopen_f')
call h5dget_space_f(dset_id, filespace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5dget_space_f')
!--------------------------------------------------------------------------------------------------
! select a hyperslab (the portion of the current process) in the file
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, localShape, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5sselect_hyperslab_f')
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, localShape, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5sselect_hyperslab_f')
end subroutine initialize_read
@ -1828,19 +1828,19 @@ end subroutine initialize_read
!--------------------------------------------------------------------------------------------------
subroutine finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
integer(HID_T), intent(in) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
integer :: hdferr
integer(HID_T), intent(in) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
integer :: hdferr
call h5pclose_f(plist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: plist_id')
call h5pclose_f(aplist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: aplist_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: h5dclose_f')
call h5sclose_f(filespace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: h5sclose_f/filespace_id')
call h5sclose_f(memspace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: h5sclose_f/memspace_id')
call h5pclose_f(plist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: plist_id')
call h5pclose_f(aplist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: aplist_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: h5dclose_f')
call h5sclose_f(filespace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: h5sclose_f/filespace_id')
call h5sclose_f(memspace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: h5sclose_f/memspace_id')
end subroutine finalize_read
@ -1852,60 +1852,60 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
myStart, totalShape, &
loc_id,myShape,datasetName,datatype,parallel)
integer(HID_T), intent(in) :: loc_id !< file or group handle
character(len=*), intent(in) :: datasetName !< name of the dataset in the file
logical, intent(in) :: parallel
integer(HID_T), intent(in) :: datatype
integer(HSIZE_T), intent(in), dimension(:) :: &
myShape
integer(HSIZE_T), intent(out), dimension(size(myShape,1)):: &
myStart, &
totalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id
integer(HID_T), intent(in) :: loc_id !< file or group handle
character(len=*), intent(in) :: datasetName !< name of the dataset in the file
logical, intent(in) :: parallel
integer(HID_T), intent(in) :: datatype
integer(HSIZE_T), intent(in), dimension(:) :: &
myShape
integer(HSIZE_T), intent(out), dimension(size(myShape,1)):: &
myStart, &
totalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id
integer, dimension(worldsize) :: &
writeSize !< contribution of all processes
integer :: ierr
integer :: hdferr
integer, dimension(worldsize) :: &
writeSize !< contribution of all processes
integer :: ierr
integer :: hdferr
!-------------------------------------------------------------------------------------------------
! creating a property list for transfer properties (is collective when reading in parallel)
call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5pcreate_f')
call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5pcreate_f')
#ifdef PETSc
if (parallel) then
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5pset_dxpl_mpio_f')
endif
if (parallel) then
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5pset_dxpl_mpio_f')
endif
#endif
!--------------------------------------------------------------------------------------------------
! determine the global data layout among all processes
writeSize = 0
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)))
writeSize = 0
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)))
#ifdef PETSc
if (parallel) then
call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if (ierr /= 0) call IO_error(894,ext_msg='initialize_write: MPI_allreduce')
endif
if (parallel) then
call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if (ierr /= 0) call IO_error(894,ext_msg='initialize_write: MPI_allreduce')
endif
#endif
myStart = int(0,HSIZE_T)
myStart(ubound(myStart)) = int(sum(writeSize(1:worldrank)),HSIZE_T)
totalShape = [myShape(1:ubound(myShape,1)-1),int(sum(writeSize),HSIZE_T)]
myStart = int(0,HSIZE_T)
myStart(ubound(myStart)) = int(sum(writeSize(1:worldrank)),HSIZE_T)
totalShape = [myShape(1:ubound(myShape,1)-1),int(sum(writeSize),HSIZE_T)]
!--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape) and in file (global shape)
call h5screate_simple_f(size(myShape), myShape, memspace_id, hdferr, myShape)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5dopen_f')
call h5screate_simple_f(size(totalShape), totalShape, filespace_id, hdferr, totalShape)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5dget_space_f')
call h5screate_simple_f(size(myShape), myShape, memspace_id, hdferr, myShape)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5dopen_f')
call h5screate_simple_f(size(totalShape), totalShape, filespace_id, hdferr, totalShape)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5dget_space_f')
!--------------------------------------------------------------------------------------------------
! create dataset in the file and select a hyperslab from it (the portion of the current process)
call h5dcreate_f(loc_id, trim(datasetName), datatype, filespace_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5dcreate_f')
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, myShape, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5sselect_hyperslab_f')
call h5dcreate_f(loc_id, trim(datasetName), datatype, filespace_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5dcreate_f')
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, myShape, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5sselect_hyperslab_f')
end subroutine initialize_write
@ -1915,19 +1915,19 @@ end subroutine initialize_write
!--------------------------------------------------------------------------------------------------
subroutine finalize_write(plist_id, dset_id, filespace_id, memspace_id)
integer(HID_T), intent(in) :: dset_id, filespace_id, memspace_id, plist_id
integer :: hdferr
integer(HID_T), intent(in) :: dset_id, filespace_id, memspace_id, plist_id
integer :: hdferr
call h5pclose_f(plist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: plist_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: h5dclose_f')
call h5sclose_f(filespace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: h5sclose_f/filespace_id')
call h5sclose_f(memspace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: h5sclose_f/memspace_id')
call h5pclose_f(plist_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: plist_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: h5dclose_f')
call h5sclose_f(filespace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: h5sclose_f/filespace_id')
call h5sclose_f(memspace_id, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: h5sclose_f/memspace_id')
end subroutine finalize_write
#endif
end module HDF5_Utilities

File diff suppressed because it is too large Load Diff

View File

@ -26,19 +26,12 @@ module config
config_numerics, &
config_debug
!ToDo: bad names (how should one know that those variables are defined in config?)
character(len=64), dimension(:), allocatable, public, protected :: &
phase_name, & !< name of each phase
homogenization_name, & !< name of each homogenization
crystallite_name, & !< name of each crystallite setting
microstructure_name, & !< name of each microstructure
texture_name !< name of each texture
! ToDo: Remove, use size(config_phase) etc
integer, public, protected :: &
material_Nphase, & !< number of phases
material_Nhomogenization !< number of homogenizations
config_name_phase, & !< name of each phase
config_name_homogenization, & !< name of each homogenization
config_name_crystallite, & !< name of each crystallite setting
config_name_microstructure, & !< name of each microstructure
config_name_texture !< name of each texture
public :: &
config_init, &
@ -81,36 +74,33 @@ subroutine config_init
select case (trim(part))
case (trim('phase'))
call parse_materialConfig(phase_name,config_phase,line,fileContent(i+1:))
call parse_materialConfig(config_name_phase,config_phase,line,fileContent(i+1:))
if (verbose) write(6,'(a)') ' Phase parsed'; flush(6)
case (trim('microstructure'))
call parse_materialConfig(microstructure_name,config_microstructure,line,fileContent(i+1:))
call parse_materialConfig(config_name_microstructure,config_microstructure,line,fileContent(i+1:))
if (verbose) write(6,'(a)') ' Microstructure parsed'; flush(6)
case (trim('crystallite'))
call parse_materialConfig(crystallite_name,config_crystallite,line,fileContent(i+1:))
call parse_materialConfig(config_name_crystallite,config_crystallite,line,fileContent(i+1:))
if (verbose) write(6,'(a)') ' Crystallite parsed'; flush(6)
case (trim('homogenization'))
call parse_materialConfig(homogenization_name,config_homogenization,line,fileContent(i+1:))
call parse_materialConfig(config_name_homogenization,config_homogenization,line,fileContent(i+1:))
if (verbose) write(6,'(a)') ' Homogenization parsed'; flush(6)
case (trim('texture'))
call parse_materialConfig(texture_name,config_texture,line,fileContent(i+1:))
call parse_materialConfig(config_name_texture,config_texture,line,fileContent(i+1:))
if (verbose) write(6,'(a)') ' Texture parsed'; flush(6)
end select
enddo
material_Nhomogenization = size(config_homogenization)
material_Nphase = size(config_phase)
if (material_Nhomogenization < 1) call IO_error(160,ext_msg='<homogenization>')
if (size(config_homogenization) < 1) call IO_error(160,ext_msg='<homogenization>')
if (size(config_microstructure) < 1) call IO_error(160,ext_msg='<microstructure>')
if (size(config_crystallite) < 1) call IO_error(160,ext_msg='<crystallite>')
if (material_Nphase < 1) call IO_error(160,ext_msg='<phase>')
if (size(config_phase) < 1) call IO_error(160,ext_msg='<phase>')
if (size(config_texture) < 1) call IO_error(160,ext_msg='<texture>')

File diff suppressed because it is too large Load Diff

View File

@ -364,7 +364,7 @@ subroutine crystallite_init
do r = 1,size(config_crystallite)
if (any(microstructure_crystallite(discretization_microstructureAt) == r)) then
write(FILEUNIT,'(/,a,/)') '['//trim(crystallite_name(r))//']'
write(FILEUNIT,'(/,a,/)') '['//trim(config_name_crystallite(r))//']'
do o = 1,crystallite_Noutput(r)
write(FILEUNIT,'(a,i4)') trim(crystallite_output(o,r))//char(9),crystallite_sizePostResult(o,r)
enddo
@ -386,7 +386,7 @@ subroutine crystallite_init
crystallite_Fp0(1:3,1:3,c,i,e) = math_EulerToR(material_EulerAngles(1:3,c,i,e)) ! plastic def gradient reflects init orientation
crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e)
crystallite_F0(1:3,1:3,c,i,e) = math_I3
crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phase(c,i,e))
crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phaseAt(c,e))
crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), &
crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration
crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e)
@ -483,12 +483,12 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then
plasticState (phaseAt(c,i,e))%subState0( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%partionedState0(:,phasememberAt(c,i,e))
plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%partionedState0(:,material_phaseMemberAt(c,i,e))
do s = 1, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(s)%subState0( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(s)%partionedState0(:,phasememberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState(material_phaseAt(c,e))%p(s)%subState0( :,material_phaseMemberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phaseMemberAt(c,i,e))
enddo
crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e)
crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_partionedLp0(1:3,1:3,c,i,e)
@ -543,11 +543,11 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e)
crystallite_subS0 (1:3,1:3,c,i,e) = crystallite_S (1:3,1:3,c,i,e)
!if abbrevation, make c and p private in omp
plasticState( phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e)) &
= plasticState(phaseAt(c,i,e))%state( :,phasememberAt(c,i,e))
do s = 1, phase_Nsources(phaseAt(c,i,e))
sourceState( phaseAt(c,i,e))%p(s)%subState0(:,phasememberAt(c,i,e)) &
= sourceState(phaseAt(c,i,e))%p(s)%state( :,phasememberAt(c,i,e))
plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) &
= plasticState(material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState( material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e)) &
= sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e))
enddo
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0 &
@ -572,11 +572,11 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_Lp (1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e)
crystallite_Li (1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e)
endif
plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) &
= plasticState(phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e))
do s = 1, phase_Nsources(phaseAt(c,i,e))
sourceState( phaseAt(c,i,e))%p(s)%state( :,phasememberAt(c,i,e)) &
= sourceState(phaseAt(c,i,e))%p(s)%subState0(:,phasememberAt(c,i,e))
plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) &
= plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState( material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e)) &
= sourceState(material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e))
enddo
! cant restore dotState here, since not yet calculated in first cutback after initialization
@ -839,7 +839,7 @@ subroutine crystallite_orientations
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
if (plasticState(material_phase(1,i,e))%nonLocal) & ! if nonlocal model
if (plasticState(material_phaseAt(1,e))%nonLocal) & ! if nonlocal model
call plastic_nonlocal_updateCompatibility(crystallite_orientation,i,e)
enddo; enddo
!$OMP END PARALLEL DO
@ -873,106 +873,106 @@ end function crystallite_push33ToRef
!--------------------------------------------------------------------------------------------------
function crystallite_postResults(ipc, ip, el)
integer, intent(in):: &
el, & !< element index
ip, & !< integration point index
ipc !< grain index
integer, intent(in):: &
el, & !< element index
ip, & !< integration point index
ipc !< grain index
real(pReal), dimension(1+crystallite_sizePostResults(microstructure_crystallite(discretization_microstructureAt(el))) + &
1+plasticState(material_phase(ipc,ip,el))%sizePostResults + &
sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: &
crystallite_postResults
integer :: &
o, &
c, &
crystID, &
mySize, &
n
type(rotation) :: rot
real(pReal), dimension(1+crystallite_sizePostResults(microstructure_crystallite(discretization_microstructureAt(el))) + &
1+plasticState(material_phaseAt(ipc,el))%sizePostResults + &
sum(sourceState(material_phaseAt(ipc,el))%p(:)%sizePostResults)) :: &
crystallite_postResults
integer :: &
o, &
c, &
crystID, &
mySize, &
n
type(rotation) :: rot
crystID = microstructure_crystallite(discretization_microstructureAt(el))
crystID = microstructure_crystallite(discretization_microstructureAt(el))
crystallite_postResults = 0.0_pReal
crystallite_postResults(1) = real(crystallite_sizePostResults(crystID),pReal) ! header-like information (length)
c = 1
crystallite_postResults = 0.0_pReal
crystallite_postResults(1) = real(crystallite_sizePostResults(crystID),pReal) ! header-like information (length)
c = 1
do o = 1,crystallite_Noutput(crystID)
mySize = 0
select case(crystallite_outputID(o,crystID))
case (phase_ID)
mySize = 1
crystallite_postResults(c+1) = real(material_phase(ipc,ip,el),pReal) ! phaseID of grain
case (texture_ID)
mySize = 1
crystallite_postResults(c+1) = real(material_texture(ipc,ip,el),pReal) ! textureID of grain
case (orientation_ID)
mySize = 4
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(ipc,ip,el)%asQuaternion()
do o = 1,crystallite_Noutput(crystID)
mySize = 0
select case(crystallite_outputID(o,crystID))
case (phase_ID)
mySize = 1
crystallite_postResults(c+1) = real(material_phaseAt(ipc,el),pReal) ! phaseID of grain
case (texture_ID)
mySize = 1
crystallite_postResults(c+1) = real(material_texture(ipc,ip,el),pReal) ! textureID of grain
case (orientation_ID)
mySize = 4
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(ipc,ip,el)%asQuaternion()
case (grainrotation_ID)
rot = crystallite_orientation0(ipc,ip,el)%misorientation(crystallite_orientation(ipc,ip,el))
mySize = 4
crystallite_postResults(c+1:c+mySize) = rot%asAxisAnglePair()
crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree
case (grainrotation_ID)
rot = crystallite_orientation0(ipc,ip,el)%misorientation(crystallite_orientation(ipc,ip,el))
mySize = 4
crystallite_postResults(c+1:c+mySize) = rot%asAxisAnglePair()
crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree
! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33
! thus row index i is slow, while column index j is fast. reminder: "row is slow"
case (defgrad_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)),[mySize])
case (fe_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)),[mySize])
case (fp_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Fp(1:3,1:3,ipc,ip,el)),[mySize])
case (fi_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Fi(1:3,1:3,ipc,ip,el)),[mySize])
case (lp_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Lp(1:3,1:3,ipc,ip,el)),[mySize])
case (li_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Li(1:3,1:3,ipc,ip,el)),[mySize])
case (p_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_P(1:3,1:3,ipc,ip,el)),[mySize])
case (s_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(crystallite_S(1:3,1:3,ipc,ip,el),[mySize])
case (elasmatrix_ID)
mySize = 36
crystallite_postResults(c+1:c+mySize) = reshape(constitutive_homogenizedC(ipc,ip,el),[mySize])
case(neighboringelement_ID)
mySize = nIPneighbors
crystallite_postResults(c+1:c+mySize) = 0.0_pReal
forall (n = 1:mySize) &
crystallite_postResults(c+n) = real(IPneighborhood(1,n,ip,el),pReal)
case(neighboringip_ID)
mySize = nIPneighbors
crystallite_postResults(c+1:c+mySize) = 0.0_pReal
forall (n = 1:mySize) &
crystallite_postResults(c+n) = real(IPneighborhood(2,n,ip,el),pReal)
end select
c = c + mySize
enddo
case (defgrad_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)),[mySize])
case (fe_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)),[mySize])
case (fp_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Fp(1:3,1:3,ipc,ip,el)),[mySize])
case (fi_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Fi(1:3,1:3,ipc,ip,el)),[mySize])
case (lp_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Lp(1:3,1:3,ipc,ip,el)),[mySize])
case (li_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Li(1:3,1:3,ipc,ip,el)),[mySize])
case (p_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_P(1:3,1:3,ipc,ip,el)),[mySize])
case (s_ID)
mySize = 9
crystallite_postResults(c+1:c+mySize) = &
reshape(crystallite_S(1:3,1:3,ipc,ip,el),[mySize])
case (elasmatrix_ID)
mySize = 36
crystallite_postResults(c+1:c+mySize) = reshape(constitutive_homogenizedC(ipc,ip,el),[mySize])
case(neighboringelement_ID)
mySize = nIPneighbors
crystallite_postResults(c+1:c+mySize) = 0.0_pReal
forall (n = 1:mySize) &
crystallite_postResults(c+n) = real(IPneighborhood(1,n,ip,el),pReal)
case(neighboringip_ID)
mySize = nIPneighbors
crystallite_postResults(c+1:c+mySize) = 0.0_pReal
forall (n = 1:mySize) &
crystallite_postResults(c+n) = real(IPneighborhood(2,n,ip,el),pReal)
end select
c = c + mySize
enddo
crystallite_postResults(c+1) = real(plasticState(material_phase(ipc,ip,el))%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)
crystallite_postResults(c+1) = real(plasticState(material_phaseAt(ipc,el))%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
@ -982,9 +982,6 @@ end function crystallite_postResults
!--------------------------------------------------------------------------------------------------
subroutine crystallite_results
#if defined(PETSc) || defined(DAMASK_HDF5)
use config, only: &
config_name_phase => phase_name ! anticipate logical name
integer :: p,o
real(pReal), allocatable, dimension(:,:,:) :: selected_tensors
type(rotation), allocatable, dimension(:) :: selected_rotations
@ -1049,13 +1046,13 @@ subroutine crystallite_results
'crystal orientation as quaternion',lattice_label)
end select
enddo
enddo
enddo
contains
contains
!--------------------------------------------------------------------------------------------------
!> @brief select tensors for output
!--------------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------------
!> @brief select tensors for output
!------------------------------------------------------------------------------------------------
function select_tensors(dataset,instance)
integer, intent(in) :: instance
@ -1094,7 +1091,7 @@ subroutine crystallite_results
j=0
do e = 1, size(material_phaseAt,2)
do i = 1, homogenization_maxNgrains !ToDo: this needs to be changed for varying Ngrains
do i = 1, homogenization_maxNgrains !ToDo: this needs to be changed for varying Ngrains
do c = 1, size(material_phaseAt,1)
if (material_phaseAt(c,e) == instance) then
j = j + 1
@ -1106,8 +1103,6 @@ subroutine crystallite_results
end function select_rotations
#endif
end subroutine crystallite_results
@ -1195,10 +1190,8 @@ logical function integrateStress(ipc,ip,el,timeFraction)
Fg_new = crystallite_subF(1:3,1:3,ipc,ip,el)
endif
!* feed local variables
Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! ... and take it as first guess
Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! ... and take it as first guess
Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess
Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess
Liguess_old = Liguess
invFp_current = math_inv33(crystallite_subFp0(1:3,1:3,ipc,ip,el))
@ -1555,7 +1548,7 @@ subroutine integrateStateFPI
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
plasticState(p)%previousDotState2(:,c) = merge(plasticState(p)%previousDotState(:,c),&
0.0_pReal,&
@ -1583,7 +1576,7 @@ subroutine integrateStateFPI
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
zeta = damper(plasticState(p)%dotState (:,c), &
@ -1746,7 +1739,7 @@ subroutine integrateStateAdaptiveEuler
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) &
@ -1775,7 +1768,7 @@ subroutine integrateStateAdaptiveEuler
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
residuum_plastic(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) &
@ -1835,7 +1828,7 @@ subroutine integrateStateRK4
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
plasticState(p)%RK4dotState(:,c) = WEIGHT(n)*plasticState(p)%dotState(:,c) &
+ merge(plasticState(p)%RK4dotState(:,c),0.0_pReal,n>1)
@ -1926,7 +1919,7 @@ subroutine integrateStateRKCK45
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e); cc = phasememberAt(g,i,e)
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc)
plasticState(p)%dotState(:,cc) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,cc)
@ -1966,7 +1959,7 @@ subroutine integrateStateRKCK45
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e); cc = phasememberAt(g,i,e)
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
@ -2005,7 +1998,7 @@ subroutine integrateStateRKCK45
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e); cc = phasememberAt(g,i,e)
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
@ -2163,7 +2156,7 @@ subroutine update_state(timeFraction)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
mySize = plasticState(p)%sizeDotState
plasticState(p)%state(1:mySize,c) = plasticState(p)%subState0(1:mySize,c) &
@ -2214,7 +2207,7 @@ subroutine update_dotState(timeFraction)
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
crystallite_subdt(g,i,e)*timeFraction, g,i,e)
p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1, phase_Nsources(p)
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
@ -2259,7 +2252,7 @@ subroutine update_deltaState
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fi(1:3,1:3,g,i,e), &
g,i,e)
p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
myOffset = plasticState(p)%offsetDeltaState
mySize = plasticState(p)%sizeDeltaState
NaN = any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))
@ -2311,8 +2304,8 @@ logical function stateJump(ipc,ip,el)
myOffset, &
mySize
c = phasememberAt(ipc,ip,el)
p = phaseAt(ipc,ip,el)
c = material_phaseMemberAt(ipc,ip,el)
p = material_phaseAt(ipc,el)
call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), &
crystallite_Fe(1:3,1:3,ipc,ip,el), &

View File

@ -178,8 +178,8 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)

View File

@ -144,8 +144,8 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
@ -194,7 +194,7 @@ function damage_nonlocal_getDiffusion33(ip,el)
damage_nonlocal_getDiffusion33 = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog)
damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + &
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phase(grain,ip,el)))
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phaseAt(grain,el)))
enddo
damage_nonlocal_getDiffusion33 = &
@ -217,7 +217,7 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
damage_nonlocal_getMobility = 0.0_pReal
do ipc = 1, homogenization_Ngrains(material_homogenizationAt(el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phaseAt(ipc,el))
enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility/&

View File

@ -6,70 +6,70 @@
!> @brief Reading in and interpretating the debugging settings for the various modules
!--------------------------------------------------------------------------------------------------
module debug
use prec
use IO
use prec
use IO
implicit none
private
implicit none
private
integer, parameter, public :: &
debug_LEVELSELECTIVE = 2**0, &
debug_LEVELBASIC = 2**1, &
debug_LEVELEXTENSIVE = 2**2
integer, parameter, private :: &
debug_MAXGENERAL = debug_LEVELEXTENSIVE ! must be set to the last bitcode used by (potentially) all debug types
integer, parameter, public :: &
debug_SPECTRALRESTART = debug_MAXGENERAL*2**1, &
debug_SPECTRALFFTW = debug_MAXGENERAL*2**2, &
debug_SPECTRALDIVERGENCE = debug_MAXGENERAL*2**3, &
debug_SPECTRALROTATION = debug_MAXGENERAL*2**4, &
debug_SPECTRALPETSC = debug_MAXGENERAL*2**5
integer, parameter, public :: &
debug_LEVELSELECTIVE = 2**0, &
debug_LEVELBASIC = 2**1, &
debug_LEVELEXTENSIVE = 2**2
integer, parameter, private :: &
debug_MAXGENERAL = debug_LEVELEXTENSIVE ! must be set to the last bitcode used by (potentially) all debug types
integer, parameter, public :: &
debug_SPECTRALRESTART = debug_MAXGENERAL*2**1, &
debug_SPECTRALFFTW = debug_MAXGENERAL*2**2, &
debug_SPECTRALDIVERGENCE = debug_MAXGENERAL*2**3, &
debug_SPECTRALROTATION = debug_MAXGENERAL*2**4, &
debug_SPECTRALPETSC = debug_MAXGENERAL*2**5
integer, parameter, public :: &
debug_DEBUG = 1, &
debug_MATH = 2, &
debug_FESOLVING = 3, &
debug_MESH = 4, & !< stores debug level for mesh part of DAMASK bitwise coded
debug_MATERIAL = 5, & !< stores debug level for material part of DAMASK bitwise coded
debug_LATTICE = 6, & !< stores debug level for lattice part of DAMASK bitwise coded
debug_CONSTITUTIVE = 7, & !< stores debug level for constitutive part of DAMASK bitwise coded
debug_CRYSTALLITE = 8, &
debug_HOMOGENIZATION = 9, &
debug_CPFEM = 10, &
debug_SPECTRAL = 11, &
debug_MARC = 12, &
debug_ABAQUS = 13
integer, parameter, private :: &
debug_MAXNTYPE = debug_ABAQUS !< must be set to the maximum defined debug type
integer, parameter, public :: &
debug_DEBUG = 1, &
debug_MATH = 2, &
debug_FESOLVING = 3, &
debug_MESH = 4, & !< stores debug level for mesh part of DAMASK bitwise coded
debug_MATERIAL = 5, & !< stores debug level for material part of DAMASK bitwise coded
debug_LATTICE = 6, & !< stores debug level for lattice part of DAMASK bitwise coded
debug_CONSTITUTIVE = 7, & !< stores debug level for constitutive part of DAMASK bitwise coded
debug_CRYSTALLITE = 8, &
debug_HOMOGENIZATION = 9, &
debug_CPFEM = 10, &
debug_SPECTRAL = 11, &
debug_MARC = 12, &
debug_ABAQUS = 13
integer, parameter, private :: &
debug_MAXNTYPE = debug_ABAQUS !< must be set to the maximum defined debug type
integer,protected, dimension(debug_maxNtype+2), public :: & ! specific ones, and 2 for "all" and "other"
debug_level = 0
integer,protected, dimension(debug_maxNtype+2), public :: & ! specific ones, and 2 for "all" and "other"
debug_level = 0
integer, protected, public :: &
debug_e = 1, &
debug_i = 1, &
debug_g = 1
integer, protected, public :: &
debug_e = 1, &
debug_i = 1, &
debug_g = 1
integer, dimension(2), public :: &
debug_stressMaxLocation = 0, &
debug_stressMinLocation = 0, &
debug_jacobianMaxLocation = 0, &
debug_jacobianMinLocation = 0
integer, dimension(2), public :: &
debug_stressMaxLocation = 0, &
debug_stressMinLocation = 0, &
debug_jacobianMaxLocation = 0, &
debug_jacobianMinLocation = 0
real(pReal), public :: &
debug_stressMax = -huge(1.0_pReal), &
debug_stressMin = huge(1.0_pReal), &
debug_jacobianMax = -huge(1.0_pReal), &
debug_jacobianMin = huge(1.0_pReal)
real(pReal), public :: &
debug_stressMax = -huge(1.0_pReal), &
debug_stressMin = huge(1.0_pReal), &
debug_jacobianMax = -huge(1.0_pReal), &
debug_jacobianMin = huge(1.0_pReal)
#ifdef PETSc
character(len=1024), parameter, public :: &
PETSCDEBUG = ' -snes_view -snes_monitor '
character(len=1024), parameter, public :: &
PETSCDEBUG = ' -snes_view -snes_monitor '
#endif
public :: debug_init, &
debug_reset, &
debug_info
public :: debug_init, &
debug_reset, &
debug_info
contains
@ -79,111 +79,111 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine debug_init
character(len=pStringLen), dimension(:), allocatable :: fileContent
character(len=pStringLen), dimension(:), allocatable :: fileContent
integer :: i, what, j
integer, allocatable, dimension(:) :: chunkPos
character(len=pStringLen) :: tag, line
logical :: fexist
integer :: i, what, j
integer, allocatable, dimension(:) :: chunkPos
character(len=pStringLen) :: tag, line
logical :: fexist
write(6,'(/,a)') ' <<<+- debug init -+>>>'
write(6,'(/,a)') ' <<<+- debug init -+>>>'
#ifdef DEBUG
write(6,'(a)') achar(27)//'[31m <<<+- DEBUG version -+>>>'//achar(27)//'[0m'
write(6,'(a)') achar(27)//'[31m <<<+- DEBUG version -+>>>'//achar(27)//'[0m'
#endif
inquire(file='debug.config', exist=fexist)
inquire(file='debug.config', exist=fexist)
fileExists: if (fexist) then
fileContent = IO_read_ASCII('debug.config')
do j=1, size(fileContent)
line = fileContent(j)
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1)) ! extract key
select case(tag)
case ('element','e','el')
debug_e = IO_intValue(line,chunkPos,2)
case ('integrationpoint','i','ip')
debug_i = IO_intValue(line,chunkPos,2)
case ('grain','g','gr')
debug_g = IO_intValue(line,chunkPos,2)
end select
fileExists: if (fexist) then
fileContent = IO_read_ASCII('debug.config')
do j=1, size(fileContent)
line = fileContent(j)
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1)) ! extract key
select case(tag)
case ('element','e','el')
debug_e = IO_intValue(line,chunkPos,2)
case ('integrationpoint','i','ip')
debug_i = IO_intValue(line,chunkPos,2)
case ('grain','g','gr')
debug_g = IO_intValue(line,chunkPos,2)
end select
what = 0
select case(tag)
case ('debug')
what = debug_DEBUG
case ('math')
what = debug_MATH
case ('fesolving', 'fe')
what = debug_FESOLVING
case ('mesh')
what = debug_MESH
case ('material')
what = debug_MATERIAL
case ('lattice')
what = debug_LATTICE
case ('constitutive')
what = debug_CONSTITUTIVE
case ('crystallite')
what = debug_CRYSTALLITE
case ('homogenization')
what = debug_HOMOGENIZATION
case ('cpfem')
what = debug_CPFEM
case ('spectral')
what = debug_SPECTRAL
case ('marc')
what = debug_MARC
case ('abaqus')
what = debug_ABAQUS
case ('all')
what = debug_MAXNTYPE + 1
case ('other')
what = debug_MAXNTYPE + 2
end select
if (what /= 0) then
do i = 2, chunkPos(1)
select case(IO_lc(IO_stringValue(line,chunkPos,i)))
case('basic')
debug_level(what) = ior(debug_level(what), debug_LEVELBASIC)
case('extensive')
debug_level(what) = ior(debug_level(what), debug_LEVELEXTENSIVE)
case('selective')
debug_level(what) = ior(debug_level(what), debug_LEVELSELECTIVE)
case('restart')
debug_level(what) = ior(debug_level(what), debug_SPECTRALRESTART)
case('fft','fftw')
debug_level(what) = ior(debug_level(what), debug_SPECTRALFFTW)
case('divergence')
debug_level(what) = ior(debug_level(what), debug_SPECTRALDIVERGENCE)
case('rotation')
debug_level(what) = ior(debug_level(what), debug_SPECTRALROTATION)
case('petsc')
debug_level(what) = ior(debug_level(what), debug_SPECTRALPETSC)
end select
enddo
endif
enddo
what = 0
select case(tag)
case ('debug')
what = debug_DEBUG
case ('math')
what = debug_MATH
case ('fesolving', 'fe')
what = debug_FESOLVING
case ('mesh')
what = debug_MESH
case ('material')
what = debug_MATERIAL
case ('lattice')
what = debug_LATTICE
case ('constitutive')
what = debug_CONSTITUTIVE
case ('crystallite')
what = debug_CRYSTALLITE
case ('homogenization')
what = debug_HOMOGENIZATION
case ('cpfem')
what = debug_CPFEM
case ('spectral')
what = debug_SPECTRAL
case ('marc')
what = debug_MARC
case ('abaqus')
what = debug_ABAQUS
case ('all')
what = debug_MAXNTYPE + 1
case ('other')
what = debug_MAXNTYPE + 2
end select
if (what /= 0) then
do i = 2, chunkPos(1)
select case(IO_lc(IO_stringValue(line,chunkPos,i)))
case('basic')
debug_level(what) = ior(debug_level(what), debug_LEVELBASIC)
case('extensive')
debug_level(what) = ior(debug_level(what), debug_LEVELEXTENSIVE)
case('selective')
debug_level(what) = ior(debug_level(what), debug_LEVELSELECTIVE)
case('restart')
debug_level(what) = ior(debug_level(what), debug_SPECTRALRESTART)
case('fft','fftw')
debug_level(what) = ior(debug_level(what), debug_SPECTRALFFTW)
case('divergence')
debug_level(what) = ior(debug_level(what), debug_SPECTRALDIVERGENCE)
case('rotation')
debug_level(what) = ior(debug_level(what), debug_SPECTRALROTATION)
case('petsc')
debug_level(what) = ior(debug_level(what), debug_SPECTRALPETSC)
end select
enddo
endif
enddo
do i = 1, debug_maxNtype
if (debug_level(i) == 0) &
debug_level(i) = ior(debug_level(i), debug_level(debug_MAXNTYPE + 2)) ! fill undefined debug types with levels specified by "other"
do i = 1, debug_maxNtype
if (debug_level(i) == 0) &
debug_level(i) = ior(debug_level(i), debug_level(debug_MAXNTYPE + 2)) ! fill undefined debug types with levels specified by "other"
debug_level(i) = ior(debug_level(i), debug_level(debug_MAXNTYPE + 1)) ! fill all debug types with levels specified by "all"
enddo
debug_level(i) = ior(debug_level(i), debug_level(debug_MAXNTYPE + 1)) ! fill all debug types with levels specified by "all"
enddo
if (iand(debug_level(debug_debug),debug_LEVELBASIC) /= 0) &
write(6,'(a,/)') ' using values from config file'
else fileExists
if (iand(debug_level(debug_debug),debug_LEVELBASIC) /= 0) &
write(6,'(a,/)') ' using standard values'
endif fileExists
if (iand(debug_level(debug_debug),debug_LEVELBASIC) /= 0) &
write(6,'(a,/)') ' using values from config file'
else fileExists
if (iand(debug_level(debug_debug),debug_LEVELBASIC) /= 0) &
write(6,'(a,/)') ' using standard values'
endif fileExists
!--------------------------------------------------------------------------------------------------
! output switched on (debug level for debug must be extensive)
if (iand(debug_level(debug_debug),debug_LEVELEXTENSIVE) /= 0) then
if (iand(debug_level(debug_debug),debug_LEVELEXTENSIVE) /= 0) then
do i = 1, debug_MAXNTYPE
select case(i)
case (debug_DEBUG)
@ -231,7 +231,7 @@ subroutine debug_init
if(iand(debug_level(i),debug_SPECTRALPETSC) /= 0) write(6,'(a)') ' PETSc'
endif
enddo
endif
endif
end subroutine debug_init
@ -241,14 +241,14 @@ end subroutine debug_init
!--------------------------------------------------------------------------------------------------
subroutine debug_reset
debug_stressMaxLocation = 0
debug_stressMinLocation = 0
debug_jacobianMaxLocation = 0
debug_jacobianMinLocation = 0
debug_stressMax = -huge(1.0_pReal)
debug_stressMin = huge(1.0_pReal)
debug_jacobianMax = -huge(1.0_pReal)
debug_jacobianMin = huge(1.0_pReal)
debug_stressMaxLocation = 0
debug_stressMinLocation = 0
debug_jacobianMaxLocation = 0
debug_jacobianMinLocation = 0
debug_stressMax = -huge(1.0_pReal)
debug_stressMin = huge(1.0_pReal)
debug_jacobianMax = -huge(1.0_pReal)
debug_jacobianMin = huge(1.0_pReal)
end subroutine debug_reset
@ -258,18 +258,18 @@ end subroutine debug_reset
!--------------------------------------------------------------------------------------------------
subroutine debug_info
!$OMP CRITICAL (write2out)
debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 &
.and. any(debug_stressMinLocation /= 0) &
.and. any(debug_stressMaxLocation /= 0) ) then
write(6,'(2/,a,/)') ' Extreme values of returned stress and Jacobian'
write(6,'(a39)') ' value el ip'
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' Jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
endif debugOutputCPFEM
!$OMP END CRITICAL (write2out)
!$OMP CRITICAL (write2out)
debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 &
.and. any(debug_stressMinLocation /= 0) &
.and. any(debug_stressMaxLocation /= 0) ) then
write(6,'(2/,a,/)') ' Extreme values of returned stress and Jacobian'
write(6,'(a39)') ' value el ip'
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' Jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
endif debugOutputCPFEM
!$OMP END CRITICAL (write2out)
end subroutine debug_info

View File

@ -1,5 +1,6 @@
!--------------------------------------------------------------------------------------------------
!> @brief spatial discretization
!> @details serves as an abstraction layer between the different solvers and DAMASK
!--------------------------------------------------------------------------------------------------
module discretization
@ -30,10 +31,12 @@ module discretization
contains
!--------------------------------------------------------------------------------------------------
!> @brief stores the relevant information in globally accesible variables
!--------------------------------------------------------------------------------------------------
subroutine discretization_init(homogenizationAt,microstructureAt,IPcoords0,NodeCoords0)
integer, dimension(:), intent(in) :: &
integer, dimension(:), intent(in) :: &
homogenizationAt, &
microstructureAt
real(pReal), dimension(:,:), intent(in) :: &
@ -57,6 +60,9 @@ subroutine discretization_init(homogenizationAt,microstructureAt,IPcoords0,NodeC
end subroutine discretization_init
!--------------------------------------------------------------------------------------------------
!> @brief write the displacements
!--------------------------------------------------------------------------------------------------
subroutine discretization_results
#if defined(PETSc) || defined(DAMASK_HDF5)
real(pReal), dimension(:,:), allocatable :: u
@ -70,6 +76,9 @@ subroutine discretization_results
end subroutine discretization_results
!--------------------------------------------------------------------------------------------------
!> @brief stores current IP coordinates
!--------------------------------------------------------------------------------------------------
subroutine discretization_setIPcoords(IPcoords)
real(pReal), dimension(:,:), intent(in) :: IPcoords
@ -78,5 +87,4 @@ subroutine discretization_setIPcoords(IPcoords)
end subroutine discretization_setIPcoords
end module discretization

View File

@ -10,7 +10,7 @@ module element
private
!---------------------------------------------------------------------------------------------------
!> Properties of a single element (the element used in the mesh)
!> Properties of a single element
!---------------------------------------------------------------------------------------------------
type, public :: tElement
integer :: &
@ -21,11 +21,9 @@ module element
Ncellnodes, &
NcellnodesPerCell, &
nIPs, &
nIPneighbors, & ! ToDo: MD: Do all IPs in one element type have the same number of neighbors?
maxNnodeAtIP
nIPneighbors
integer, dimension(:,:), allocatable :: &
Cell, & !< intra-element (cell) nodes that constitute a cell
NnodeAtIP, &
IPneighbor, &
cellFace
integer, dimension(:,:), allocatable :: &
@ -139,21 +137,6 @@ module element
4 & ! 3D 8node
] !< number of cell nodes in a specific cell type
!integer, dimension(maxval(geomType)), parameter, private :: maxNodeAtIP = & ! Intel 16.0 complains
integer, dimension(10), parameter, private :: maxNnodeAtIP = &
[ &
3, &
1, &
1, &
2, &
4, &
1, &
1, &
8, &
1, &
4 &
] !< maximum number of parent nodes that belong to an IP for a specific type of element
!integer, dimension(maxval(CELLTYPE)), parameter, private :: NCELLNODEPERCELL = & ! Intel 16.0 complains
integer, dimension(4), parameter, private :: NCELLNODEPERCELL = &
[ &
@ -163,114 +146,6 @@ module element
8 & ! 3D 8node
] !< number of cell nodes in a specific cell type
! --------------------------------------------------------------------------------------------------
! MD: probably not needed START
integer, dimension(maxNnodeAtIP(1),nIP(1)), parameter, private :: NnodeAtIP1 = &
reshape([&
1,2,3 &
],[maxNnodeAtIP(1),nIP(1)])
integer, dimension(maxNnodeAtIP(2),nIP(2)), parameter, private :: NnodeAtIP2 = &
reshape([&
1, &
2, &
3 &
],[maxNnodeAtIP(2),nIP(2)])
integer, dimension(maxNnodeAtIP(3),nIP(3)), parameter, private :: NnodeAtIP3 = &
reshape([&
1, &
2, &
4, &
3 &
],[maxNnodeAtIP(3),nIP(3)])
integer, dimension(maxNnodeAtIP(4),nIP(4)), parameter, private :: NnodeAtIP4 = &
reshape([&
1,0, &
1,2, &
2,0, &
1,4, &
0,0, &
2,3, &
4,0, &
3,4, &
3,0 &
],[maxNnodeAtIP(4),nIP(4)])
integer, dimension(maxNnodeAtIP(5),nIP(5)), parameter, private :: NnodeAtIP5 = &
reshape([&
1,2,3,4 &
],[maxNnodeAtIP(5),nIP(5)])
integer, dimension(maxNnodeAtIP(6),nIP(6)), parameter, private :: NnodeAtIP6 = &
reshape([&
1, &
2, &
3, &
4 &
],[maxNnodeAtIP(6),nIP(6)])
integer, dimension(maxNnodeAtIP(7),nIP(7)), parameter, private :: NnodeAtIP7 = &
reshape([&
1, &
2, &
3, &
4, &
5, &
6 &
],[maxNnodeAtIP(7),nIP(7)])
integer, dimension(maxNnodeAtIP(8),nIP(8)), parameter, private :: NnodeAtIP8 = &
reshape([&
1,2,3,4,5,6,7,8 &
],[maxNnodeAtIP(8),nIP(8)])
integer, dimension(maxNnodeAtIP(9),nIP(9)), parameter, private :: NnodeAtIP9 = &
reshape([&
1, &
2, &
4, &
3, &
5, &
6, &
8, &
7 &
],[maxNnodeAtIP(9),nIP(9)])
integer, dimension(maxNnodeAtIP(10),nIP(10)), parameter, private :: NnodeAtIP10 = &
reshape([&
1,0, 0,0, &
1,2, 0,0, &
2,0, 0,0, &
1,4, 0,0, &
1,3, 2,4, &
2,3, 0,0, &
4,0, 0,0, &
3,4, 0,0, &
3,0, 0,0, &
1,5, 0,0, &
1,6, 2,5, &
2,6, 0,0, &
1,8, 4,5, &
0,0, 0,0, &
2,7, 3,6, &
4,8, 0,0, &
3,8, 4,7, &
3,7, 0,0, &
5,0, 0,0, &
5,6, 0,0, &
6,0, 0,0, &
5,8, 0,0, &
5,7, 6,8, &
6,7, 0,0, &
8,0, 0,0, &
7,8, 0,0, &
7,0, 0,0 &
],[maxNnodeAtIP(10),nIP(10)])
! *** FE_ipNeighbor ***
! is a list of the neighborhood of each IP.
! It is sorted in (local) +x,-x, +y,-y, +z,-z direction.
@ -386,15 +261,15 @@ module element
real(pReal), dimension(nNode(1),NcellNode(geomType(1))), parameter :: cellNodeParentNodeWeights1 = &
reshape(real([&
integer, dimension(nNode(1),NcellNode(geomType(1))), parameter :: cellNodeParentNodeWeights1 = &
reshape([&
1, 0, 0, &
0, 1, 0, &
0, 0, 1 &
],pReal),[nNode(1),NcellNode(geomType(1))]) ! 2D 3node 1ip
],[nNode(1),NcellNode(geomType(1))]) !< 2D 3node 1ip
real(pReal), dimension(nNode(2),NcellNode(geomType(2))), parameter :: cellNodeParentNodeWeights2 = &
reshape(real([&
integer, dimension(nNode(2),NcellNode(geomType(2))), parameter :: cellNodeParentNodeWeights2 = &
reshape([&
1, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, &
0, 0, 1, 0, 0, 0, &
@ -402,10 +277,10 @@ module element
0, 0, 0, 0, 1, 0, &
0, 0, 0, 0, 0, 1, &
1, 1, 1, 2, 2, 2 &
],pReal),[nNode(2),NcellNode(geomType(2))]) ! 2D 6node 3ip
],[nNode(2),NcellNode(geomType(2))]) !< 2D 6node 3ip
real(pReal), dimension(nNode(3),NcellNode(geomType(3))), parameter :: cellNodeParentNodeWeights3 = &
reshape(real([&
integer, dimension(nNode(3),NcellNode(geomType(3))), parameter :: cellNodeParentNodeWeights3 = &
reshape([&
1, 0, 0, 0, &
0, 1, 0, 0, &
0, 0, 1, 0, &
@ -415,10 +290,10 @@ module element
0, 0, 1, 1, &
1, 0, 0, 1, &
1, 1, 1, 1 &
],pReal),[nNode(3),NcellNode(geomType(3))]) ! 2D 6node 3ip
],[nNode(3),NcellNode(geomType(3))]) !< 2D 6node 3ip
real(pReal), dimension(nNode(4),NcellNode(geomType(4))), parameter :: cellNodeParentNodeWeights4 = &
reshape(real([&
integer, dimension(nNode(4),NcellNode(geomType(4))), parameter :: cellNodeParentNodeWeights4 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, 0, 0, &
0, 0, 1, 0, 0, 0, 0, 0, &
@ -435,10 +310,10 @@ module element
1, 4, 1, 1, 8, 8, 2, 2, &
1, 1, 4, 1, 2, 8, 8, 2, &
1, 1, 1, 4, 2, 2, 8, 8 &
],pReal),[nNode(4),NcellNode(geomType(4))]) ! 2D 8node 9ip
],[nNode(4),NcellNode(geomType(4))]) !< 2D 8node 9ip
real(pReal), dimension(nNode(5),NcellNode(geomType(5))), parameter :: cellNodeParentNodeWeights5 = &
reshape(real([&
integer, dimension(nNode(5),NcellNode(geomType(5))), parameter :: cellNodeParentNodeWeights5 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, 0, 0, &
0, 0, 1, 0, 0, 0, 0, 0, &
@ -448,18 +323,18 @@ module element
0, 0, 0, 0, 0, 0, 1, 0, &
0, 0, 0, 0, 0, 0, 0, 1, &
1, 1, 1, 1, 2, 2, 2, 2 &
],pReal),[nNode(5),NcellNode(geomType(5))]) ! 2D 8node 4ip
],[nNode(5),NcellNode(geomType(5))]) !< 2D 8node 4ip
real(pReal), dimension(nNode(6),NcellNode(geomType(6))), parameter :: cellNodeParentNodeWeights6 = &
reshape(real([&
integer, dimension(nNode(6),NcellNode(geomType(6))), parameter :: cellNodeParentNodeWeights6 = &
reshape([&
1, 0, 0, 0, &
0, 1, 0, 0, &
0, 0, 1, 0, &
0, 0, 0, 1 &
],pReal),[nNode(6),NcellNode(geomType(6))]) ! 3D 4node 1ip
],[nNode(6),NcellNode(geomType(6))]) !< 3D 4node 1ip
real(pReal), dimension(nNode(7),NcellNode(geomType(7))), parameter :: cellNodeParentNodeWeights7 = &
reshape(real([&
integer, dimension(nNode(7),NcellNode(geomType(7))), parameter :: cellNodeParentNodeWeights7 = &
reshape([&
1, 0, 0, 0, 0, &
0, 1, 0, 0, 0, &
0, 0, 1, 0, 0, &
@ -475,10 +350,10 @@ module element
0, 1, 1, 1, 0, &
1, 0, 1, 1, 0, &
0, 0, 0, 0, 1 &
],pReal),[nNode(7),NcellNode(geomType(7))]) ! 3D 5node 4ip
],[nNode(7),NcellNode(geomType(7))]) !< 3D 5node 4ip
real(pReal), dimension(nNode(8),NcellNode(geomType(8))), parameter :: cellNodeParentNodeWeights8 = &
reshape(real([&
integer, dimension(nNode(8),NcellNode(geomType(8))), parameter :: cellNodeParentNodeWeights8 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, &
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, &
@ -494,10 +369,10 @@ module element
0, 1, 1, 1, 0, 2, 0, 0, 2, 2, &
1, 0, 1, 1, 0, 0, 2, 2, 0, 2, &
3, 3, 3, 3, 4, 4, 4, 4, 4, 4 &
],pReal),[nNode(8),NcellNode(geomType(8))]) ! 3D 10node 4ip
],[nNode(8),NcellNode(geomType(8))]) !< 3D 10node 4ip
real(pReal), dimension(nNode(9),NcellNode(geomType(9))), parameter :: cellNodeParentNodeWeights9 = &
reshape(real([&
integer, dimension(nNode(9),NcellNode(geomType(9))), parameter :: cellNodeParentNodeWeights9 = &
reshape([&
1, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, &
0, 0, 1, 0, 0, 0, &
@ -519,10 +394,10 @@ module element
1, 0, 1, 1, 0, 1, &
0, 0, 0, 1, 1, 1, &
1, 1, 1, 1, 1, 1 &
],pReal),[nNode(9),NcellNode(geomType(9))]) ! 3D 6node 6ip
],[nNode(9),NcellNode(geomType(9))]) !< 3D 6node 6ip
real(pReal), dimension(nNode(10),NcellNode(geomType(10))), parameter :: cellNodeParentNodeWeights10 = &
reshape(real([&
integer, dimension(nNode(10),NcellNode(geomType(10))), parameter :: cellNodeParentNodeWeights10 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, 0, 0, &
0, 0, 1, 0, 0, 0, 0, 0, &
@ -531,10 +406,10 @@ module element
0, 0, 0, 0, 0, 1, 0, 0, &
0, 0, 0, 0, 0, 0, 1, 0, &
0, 0, 0, 0, 0, 0, 0, 1 &
],pReal),[nNode(10),NcellNode(geomType(10))]) ! 3D 8node 1ip
],[nNode(10),NcellNode(geomType(10))]) !< 3D 8node 1ip
real(pReal), dimension(nNode(11),NcellNode(geomType(11))), parameter :: cellNodeParentNodeWeights11 = &
reshape(real([&
integer, dimension(nNode(11),NcellNode(geomType(11))), parameter :: cellNodeParentNodeWeights11 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, & !
0, 1, 0, 0, 0, 0, 0, 0, & !
0, 0, 1, 0, 0, 0, 0, 0, & !
@ -562,10 +437,10 @@ module element
1, 0, 0, 1, 1, 0, 0, 1, & ! 25
0, 0, 0, 0, 1, 1, 1, 1, & !
1, 1, 1, 1, 1, 1, 1, 1 & !
],pReal),[nNode(11),NcellNode(geomType(11))]) ! 3D 8node 8ip
],[nNode(11),NcellNode(geomType(11))]) !< 3D 8node 8ip
real(pReal), dimension(nNode(12),NcellNode(geomType(12))), parameter :: cellNodeParentNodeWeights12 = &
reshape(real([&
integer, dimension(nNode(12),NcellNode(geomType(12))), parameter :: cellNodeParentNodeWeights12 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
@ -593,10 +468,10 @@ module element
1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, 2, & ! 25
0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, & !
3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 & !
],pReal),[nNode(12),NcellNode(geomType(12))]) ! 3D 20node 8ip
],[nNode(12),NcellNode(geomType(12))]) !< 3D 20node 8ip
real(pReal), dimension(nNode(13),NcellNode(geomType(13))), parameter :: cellNodeParentNodeWeights13 = &
reshape(real([&
integer, dimension(nNode(13),NcellNode(geomType(13))), parameter :: cellNodeParentNodeWeights13 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
@ -661,7 +536,7 @@ module element
4, 8, 4, 3, 8,24, 8, 4, 12,12, 4, 4, 32,32,12,12, 12,32,12, 4, & !
3, 4, 8, 4, 4, 8,24, 8, 4,12,12, 4, 12,32,32,12, 4,12,32,12, & !
4, 3, 4, 8, 8, 4, 8,24, 4, 4,12,12, 12,12,32,32, 12, 4,12,32 & !
],pReal),[nNode(13),NcellNode(geomType(13))]) ! 3D 20node 27ip
],[nNode(13),NcellNode(geomType(13))]) !< 3D 20node 27ip
integer, dimension(NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)), parameter :: CELL1 = &
@ -803,9 +678,9 @@ module element
],[NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)]) !< 3D 8node, VTK_HEXAHEDRON (12)
contains
contains
subroutine tElement_init(self,elemType)
subroutine tElement_init(self,elemType)
class(tElement) :: self
integer, intent(in) :: elemType
@ -846,50 +721,38 @@ module element
self%NcellNodes = NcellNode (self%geomType)
self%maxNnodeAtIP = maxNnodeAtIP (self%geomType)
self%nIPs = nIP (self%geomType)
self%cellType = cellType (self%geomType)
select case (self%geomType)
case(1)
self%NnodeAtIP = NnodeAtIP1
self%IPneighbor = IPneighbor1
self%cell = CELL1
case(2)
self%NnodeAtIP = NnodeAtIP2
self%IPneighbor = IPneighbor2
self%cell = CELL2
case(3)
self%NnodeAtIP = NnodeAtIP3
self%IPneighbor = IPneighbor3
self%cell = CELL3
case(4)
self%NnodeAtIP = NnodeAtIP4
self%IPneighbor = IPneighbor4
self%cell = CELL4
case(5)
self%NnodeAtIP = NnodeAtIP5
self%IPneighbor = IPneighbor5
self%cell = CELL5
case(6)
self%NnodeAtIP = NnodeAtIP6
self%IPneighbor = IPneighbor6
self%cell = CELL6
case(7)
self%NnodeAtIP = NnodeAtIP7
self%IPneighbor = IPneighbor7
self%cell = CELL7
case(8)
self%NnodeAtIP = NnodeAtIP8
self%IPneighbor = IPneighbor8
self%cell = CELL8
case(9)
self%NnodeAtIP = NnodeAtIP9
self%IPneighbor = IPneighbor9
self%cell = CELL9
case(10)
self%NnodeAtIP = NnodeAtIP10
self%IPneighbor = IPneighbor10
self%cell = CELL10
end select
@ -911,16 +774,15 @@ module element
write(6,'(/,a)') ' <<<+- element_init -+>>>'
write(6,*)' element type: ',self%elemType
write(6,*)' geom type: ',self%geomType
write(6,*)' cell type: ',self%cellType
write(6,*)' # node: ',self%Nnodes
write(6,*)' # IP: ',self%nIPs
write(6,*)' # cellnode: ',self%Ncellnodes
write(6,*)' # cellnode/cell: ',self%NcellnodesPerCell
write(6,*)' # IP neighbor: ',self%nIPneighbors
write(6,*)' max # node at IP: ',self%maxNnodeAtIP
write(6,*) ' element type: ',self%elemType
write(6,*) ' geom type: ',self%geomType
write(6,*) ' cell type: ',self%cellType
write(6,*) ' # node: ',self%Nnodes
write(6,*) ' # IP: ',self%nIPs
write(6,*) ' # cellnode: ',self%Ncellnodes
write(6,*) ' # cellnode/cell: ',self%NcellnodesPerCell
write(6,*) ' # IP neighbor: ',self%nIPneighbors
end subroutine tElement_init
end subroutine tElement_init
end module element

View File

@ -32,6 +32,7 @@ function findloc(a,v)
end function findloc
#endif
#if defined(__PGI)
!--------------------------------------------------------------------------------------------------
!> @brief substitute for the norm2 intrinsic (only for real, dimension(3) at the moment)

View File

@ -217,8 +217,7 @@ program DAMASK_spectral
case('freq','frequency','outputfreq') ! frequency of result writings
newLoadCase%outputfrequency = IO_intValue(line,chunkPos,i+1)
case('r','restart','restartwrite') ! frequency of writing restart information
newLoadCase%restartfrequency = &
max(0,IO_intValue(line,chunkPos,i+1))
newLoadCase%restartfrequency = IO_intValue(line,chunkPos,i+1)
case('guessreset','dropguessing')
newLoadCase%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('euler') ! rotation of load case given in euler angles
@ -300,7 +299,9 @@ program DAMASK_spectral
write(6,'(2x,a,i5)') 'increments: ', newLoadCase%incs
if (newLoadCase%outputfrequency < 1) errorID = 836 ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', newLoadCase%outputfrequency
write(6,'(2x,a,i5)') 'restart frequency: ', newLoadCase%restartfrequency
if (newLoadCase%restartfrequency < 1) errorID = 839 ! non-positive restart frequency
if (newLoadCase%restartfrequency < huge(0)) &
write(6,'(2x,a,i5)') 'restart frequency: ', newLoadCase%restartfrequency
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
endif reportAndCheck
loadCases = [loadCases,newLoadCase] ! load case is ok, append it
@ -336,20 +337,20 @@ program DAMASK_spectral
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) '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) '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:', restartInc ! start with writing out the previous inc
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
close(fileUnit) ! end of header
open(newunit=statUnit,file=trim(getSolverJobName())//&
'.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
@ -425,7 +426,7 @@ program DAMASK_spectral
endif
timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
skipping: if (totalIncsCounter <= restartInc) then ! not yet at restart inc?
skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc?
time = time + timeinc ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping
@ -561,8 +562,7 @@ program DAMASK_spectral
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
call CPFEM_results(totalIncsCounter,time)
endif
if ( loadCases(currentLoadCase)%restartFrequency > 0 & ! writing of restart info requested ...
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then ! ... and at frequency of writing restart information
if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then ! at frequency of writing restart information
restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif

View File

@ -181,8 +181,9 @@ subroutine grid_mech_FEM_init
!--------------------------------------------------------------------------------------------------
! init fields
restart: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file'
restartRead: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
'reading values of increment ', interface_restartInc, ' from file'
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
@ -195,10 +196,10 @@ subroutine grid_mech_FEM_init
call HDF5_read(fileHandle,u_current, 'u')
call HDF5_read(fileHandle,u_lastInc, 'u_lastInc')
elseif (restartInc == 0) then restart
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restart
endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateIPcoords(F)
call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
@ -210,12 +211,13 @@ subroutine grid_mech_FEM_init
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
restartRead: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file'
restartRead2: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
'reading more values of increment ', interface_restartInc, ' from file'
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
endif restartRead
endif restartRead2
end subroutine grid_mech_FEM_init

View File

@ -151,8 +151,9 @@ subroutine grid_mech_spectral_basic_init
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
restart: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file'
restartRead: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
' reading values of increment ', interface_restartInc, ' from file'
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
@ -163,10 +164,10 @@ subroutine grid_mech_spectral_basic_init
call HDF5_read(fileHandle,F, 'F')
call HDF5_read(fileHandle,F_lastInc, 'F_lastInc')
elseif (restartInc == 0) then restart
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restart
endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
@ -176,15 +177,16 @@ subroutine grid_mech_spectral_basic_init
math_I3) ! no rotation of boundary condition
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file'
restartRead2: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
'reading more values of increment ', interface_restartInc, ' from file'
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
fileUnit = IO_open_jobFile_binary('C_ref')
read(fileUnit) C_minMaxAvg; close(fileUnit)
endif restartRead
endif restartRead2
call utilities_updateGamma(C_minMaxAvg,.true.)

View File

@ -160,8 +160,9 @@ subroutine grid_mech_spectral_polarisation_init
F => FandF_tau( 0: 8,:,:,:)
F_tau => FandF_tau( 9:17,:,:,:)
restart: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file'
restartRead: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
' reading values of increment ', interface_restartInc, ' from file'
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
@ -174,12 +175,12 @@ subroutine grid_mech_spectral_polarisation_init
call HDF5_read(fileHandle,F_tau, 'F_tau')
call HDF5_read(fileHandle,F_tau_lastInc,'F_tau_lastInc')
elseif (restartInc == 0) then restart
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
F_tau = 2.0_pReal*F
F_tau_lastInc = 2.0_pReal*F_lastInc
endif restart
endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
@ -189,15 +190,16 @@ subroutine grid_mech_spectral_polarisation_init
math_I3) ! no rotation of boundary condition
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file'
restartRead2: if (interface_restartInc > 0) then
write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
' reading more values of increment ', interface_restartInc, ' from file'
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
fileUnit = IO_open_jobFile_binary('C_ref')
read(fileUnit) C_minMaxAvg; close(fileUnit)
endif restartRead
endif restartRead2
call utilities_updateGamma(C_minMaxAvg,.true.)
C_scale = C_minMaxAvg

View File

@ -98,7 +98,7 @@ module spectral_utilities
real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments
outputfrequency = 1, & !< frequency of result writes
restartfrequency = 0, & !< frequency of restart writes
restartfrequency = huge(0), & !< frequency of restart writes
logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -102,15 +102,15 @@ subroutine kinematics_cleavage_opening_init
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,p),& ! limit active cleavage systems per family to min of available and requested
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance))
kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether
if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether
if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
enddo
end subroutine kinematics_cleavage_opening_init
@ -138,7 +138,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
phase = material_phase(ipc,ip,el)
phase = material_phaseAt(ipc,el)
instance = kinematics_cleavage_opening_instance(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)

View File

@ -5,40 +5,40 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_slipplane_opening
use prec
use config
use IO
use debug
use math
use lattice
use material
use prec
use config
use IO
use debug
use math
use lattice
use material
implicit none
private
implicit none
private
integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance
integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance
type :: tParameters !< container type for internal constitutive parameters
integer :: &
totalNslip
integer, dimension(:), allocatable :: &
Nslip !< active number of slip systems per family
real(pReal) :: &
sdot0, &
n
real(pReal), dimension(:), allocatable :: &
critLoad
real(pReal), dimension(:,:), allocatable :: &
slip_direction, &
slip_normal, &
slip_transverse
end type tParameters
type :: tParameters !< container type for internal constitutive parameters
integer :: &
totalNslip
integer, dimension(:), allocatable :: &
Nslip !< active number of slip systems per family
real(pReal) :: &
sdot0, &
n
real(pReal), dimension(:), allocatable :: &
critLoad
real(pReal), dimension(:,:), allocatable :: &
slip_direction, &
slip_normal, &
slip_transverse
end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
kinematics_slipplane_opening_init, &
kinematics_slipplane_opening_LiAndItsTangent
public :: &
kinematics_slipplane_opening_init, &
kinematics_slipplane_opening_LiAndItsTangent
contains
@ -49,53 +49,53 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_init
integer :: maxNinstance,p,instance
integer :: maxNinstance,p,instance
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'
maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID)
if (maxNinstance == 0) return
maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID)
if (maxNinstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0)
do p = 1, size(config_phase)
kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct?
enddo
allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0)
do p = 1, size(config_phase)
kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct?
enddo
allocate(param(maxNinstance))
allocate(param(maxNinstance))
do p = 1, size(config_phase)
if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle
associate(prm => param(kinematics_slipplane_opening_instance(p)), &
config => config_phase(p))
instance = kinematics_slipplane_opening_instance(p)
prm%sdot0 = config_phase(p)%getFloat('anisoductile_sdot0')
prm%n = config_phase(p)%getFloat('anisoductile_ratesensitivity')
do p = 1, size(config_phase)
if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle
associate(prm => param(kinematics_slipplane_opening_instance(p)), &
config => config_phase(p))
instance = kinematics_slipplane_opening_instance(p)
prm%sdot0 = config_phase(p)%getFloat('anisoductile_sdot0')
prm%n = config_phase(p)%getFloat('anisoductile_ratesensitivity')
prm%Nslip = config%getInts('nslip')
prm%Nslip = config%getInts('nslip')
prm%critLoad = config_phase(p)%getFloats('anisoductile_criticalload',requiredSize=size(prm%Nslip ))
prm%critLoad = config_phase(p)%getFloats('anisoductile_criticalload',requiredSize=size(prm%Nslip ))
prm%critLoad = math_expand(prm%critLoad, prm%Nslip)
prm%critLoad = math_expand(prm%critLoad, prm%Nslip)
prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
! if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) &
! call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) &
! call IO_error(211,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) &
! call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) &
! call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) &
! call IO_error(211,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) &
! call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')')
end associate
enddo
end associate
enddo
end subroutine kinematics_slipplane_opening_init
@ -104,84 +104,84 @@ end subroutine kinematics_slipplane_opening_init
!--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
real(pReal), dimension(3,3) :: &
projection_d, projection_t, projection_n !< projection modes 3x3 tensor
integer :: &
instance, phase, &
homog, damageOffset, &
i, k, l, m, n
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
real(pReal), dimension(3,3) :: &
projection_d, projection_t, projection_n !< projection modes 3x3 tensor
integer :: &
instance, phase, &
homog, damageOffset, &
i, k, l, m, n
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
phase = material_phase(ipc,ip,el)
instance = kinematics_slipplane_opening_instance(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
phase = material_phaseAt(ipc,el)
instance = kinematics_slipplane_opening_instance(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
associate(prm => param(instance))
Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal
do i = 1, prm%totalNslip
associate(prm => param(instance))
Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal
do i = 1, prm%totalNslip
projection_d = math_outer(prm%slip_direction(1:3,i),prm%slip_normal(1:3,i))
projection_t = math_outer(prm%slip_transverse(1:3,i),prm%slip_normal(1:3,i))
projection_n = math_outer(prm%slip_normal(1:3,i),prm%slip_normal(1:3,i))
projection_d = math_outer(prm%slip_direction(1:3,i),prm%slip_normal(1:3,i))
projection_t = math_outer(prm%slip_transverse(1:3,i),prm%slip_normal(1:3,i))
projection_n = math_outer(prm%slip_normal(1:3,i),prm%slip_normal(1:3,i))
traction_d = math_mul33xx33(S,projection_d)
traction_t = math_mul33xx33(S,projection_t)
traction_n = math_mul33xx33(S,projection_n)
traction_d = math_mul33xx33(S,projection_d)
traction_t = math_mul33xx33(S,projection_t)
traction_n = math_mul33xx33(S,projection_n)
traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage
traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage
udotd = sign(1.0_pReal,traction_d)* &
prm%sdot0* &
(abs(traction_d)/traction_crit - &
abs(traction_d)/prm%critLoad(i))**prm%n
if (abs(udotd) > tol_math_check) then
Ld = Ld + udotd*projection_d
dudotd_dt = udotd*prm%n/traction_d
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotd_dt*projection_d(k,l)*projection_d(m,n)
endif
udotd = sign(1.0_pReal,traction_d)* &
prm%sdot0* &
(abs(traction_d)/traction_crit - &
abs(traction_d)/prm%critLoad(i))**prm%n
if (abs(udotd) > tol_math_check) then
Ld = Ld + udotd*projection_d
dudotd_dt = udotd*prm%n/traction_d
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotd_dt*projection_d(k,l)*projection_d(m,n)
endif
udott = sign(1.0_pReal,traction_t)* &
prm%sdot0* &
(abs(traction_t)/traction_crit - &
abs(traction_t)/prm%critLoad(i))**prm%n
if (abs(udott) > tol_math_check) then
Ld = Ld + udott*projection_t
dudott_dt = udott*prm%n/traction_t
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudott_dt*projection_t(k,l)*projection_t(m,n)
endif
udott = sign(1.0_pReal,traction_t)* &
prm%sdot0* &
(abs(traction_t)/traction_crit - &
abs(traction_t)/prm%critLoad(i))**prm%n
if (abs(udott) > tol_math_check) then
Ld = Ld + udott*projection_t
dudott_dt = udott*prm%n/traction_t
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudott_dt*projection_t(k,l)*projection_t(m,n)
endif
udotn = &
prm%sdot0* &
(max(0.0_pReal,traction_n)/traction_crit - &
max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n
if (abs(udotn) > tol_math_check) then
Ld = Ld + udotn*projection_n
dudotn_dt = udotn*prm%n/traction_n
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotn_dt*projection_n(k,l)*projection_n(m,n)
endif
enddo
udotn = &
prm%sdot0* &
(max(0.0_pReal,traction_n)/traction_crit - &
max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n
if (abs(udotn) > tol_math_check) then
Ld = Ld + udotn*projection_n
dudotn_dt = udotn*prm%n/traction_n
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotn_dt*projection_n(k,l)*projection_n(m,n)
endif
enddo
end associate
end associate
end subroutine kinematics_slipplane_opening_LiAndItsTangent

View File

@ -84,11 +84,11 @@ pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset)
kinematics_thermal_expansion_initialStrain = &
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * &
lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient
lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * &
lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient
lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * &
lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient
lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient
end function kinematics_thermal_expansion_initialStrain
@ -99,20 +99,20 @@ end function kinematics_thermal_expansion_initialStrain
subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
integer :: &
phase, &
homog, offset
real(pReal) :: &
T, TRef, TDot
phase = material_phase(ipc,ip,el)
phase = material_phaseAt(ipc,el)
homog = material_homogenizationAt(el)
offset = thermalMapping(homog)%p(ip,el)
T = temperature(homog)%p(offset)
@ -120,9 +120,9 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip,
TRef = lattice_referenceTemperature(phase)
Li = TDot * ( &
lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient
lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient
) / &
(1.0_pReal &
+ lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. &

View File

@ -98,6 +98,10 @@ module material
integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: &
damage_type !< nonlocal damage model
integer, public, protected :: &
material_Nphase, & !< number of phases
material_Nhomogenization !< number of homogenizations
integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable, public, protected :: &
phase_source, & !< active sources mechanisms of each phase
phase_kinematics, & !< active kinematic mechanisms of each phase
@ -139,10 +143,6 @@ module material
material_phaseMemberAt !< position of the element within its phase instance
! END NEW MAPPINGS
! DEPRECATED: use material_phaseAt
integer, dimension(:,:,:), allocatable, public :: &
material_phase !< phase (index) of each grain,IP,element
type(tPlasticState), allocatable, dimension(:), public :: &
plasticState
type(tSourceState), allocatable, dimension(:), public :: &
@ -180,9 +180,6 @@ module material
homogenization_active
! BEGIN DEPRECATED
integer, dimension(:,:,:), allocatable, public :: phaseAt !< phase ID of every material point (ipc,ip,el)
integer, dimension(:,:,:), allocatable, public :: phasememberAt !< memberID of given phase at every material point (ipc,ip,el)
integer, dimension(:,:,:), allocatable, public, target :: mappingHomogenization !< mapping from material points to offset in heterogenous state/field
integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
! END DEPRECATED
@ -233,25 +230,18 @@ module material
material_parseMicrostructure, &
material_parseCrystallite, &
material_parsePhase, &
material_parseTexture, &
material_populateGrains
material_parseTexture
contains
!--------------------------------------------------------------------------------------------------
!> @brief parses material configuration file
!> @details figures out if solverJobName.materialConfig is present, if not looks for
!> material.config
!--------------------------------------------------------------------------------------------------
subroutine material_init
integer, parameter :: FILEUNIT = 210
integer :: m,c,h, myDebug, myPhase, myHomog
integer :: &
g, & !< grain number
i, & !< integration point number
e !< element number
integer :: i,e,m,c,h, myDebug, myPhase, myHomog, myMicro
integer, dimension(:), allocatable :: &
CounterPhase, &
CounterHomogenization
@ -275,23 +265,27 @@ subroutine material_init
call material_parseTexture()
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Texture parsed'; flush(6)
allocate(plasticState (size(config_phase)))
allocate(sourceState (size(config_phase)))
do myPhase = 1,size(config_phase)
material_Nphase = size(config_phase)
material_Nhomogenization = size(config_homogenization)
allocate(plasticState(material_Nphase))
allocate(sourceState (material_Nphase))
do myPhase = 1,material_Nphase
allocate(sourceState(myPhase)%p(phase_Nsources(myPhase)))
enddo
allocate(homogState (size(config_homogenization)))
allocate(thermalState (size(config_homogenization)))
allocate(damageState (size(config_homogenization)))
allocate(homogState (material_Nhomogenization))
allocate(thermalState (material_Nhomogenization))
allocate(damageState (material_Nhomogenization))
allocate(thermalMapping (size(config_homogenization)))
allocate(damageMapping (size(config_homogenization)))
allocate(thermalMapping (material_Nhomogenization))
allocate(damageMapping (material_Nhomogenization))
allocate(temperature (size(config_homogenization)))
allocate(damage (size(config_homogenization)))
allocate(temperature (material_Nhomogenization))
allocate(damage (material_Nhomogenization))
allocate(temperatureRate (size(config_homogenization)))
allocate(temperatureRate (material_Nhomogenization))
do m = 1,size(config_microstructure)
if(microstructure_crystallite(m) < 1 .or. &
@ -311,17 +305,17 @@ subroutine material_init
write(6,'(/,a,/)') ' MATERIAL configuration'
write(6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains'
do h = 1,size(config_homogenization)
write(6,'(1x,a32,1x,a16,1x,i6)') homogenization_name(h),homogenization_type(h),homogenization_Ngrains(h)
write(6,'(1x,a32,1x,a16,1x,i6)') config_name_homogenization(h),homogenization_type(h),homogenization_Ngrains(h)
enddo
write(6,'(/,a14,18x,1x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents'
do m = 1,size(config_microstructure)
write(6,'(1x,a32,1x,i11,1x,i12)') microstructure_name(m), &
write(6,'(1x,a32,1x,i11,1x,i12)') config_name_microstructure(m), &
microstructure_crystallite(m), &
microstructure_Nconstituents(m)
if (microstructure_Nconstituents(m) > 0) then
do c = 1,microstructure_Nconstituents(m)
write(6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(c,m)),&
texture_name(microstructure_texture(c,m)),&
write(6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',config_name_phase(microstructure_phase(c,m)),&
config_name_texture(microstructure_texture(c,m)),&
microstructure_fraction(c,m)
enddo
write(6,*)
@ -329,10 +323,27 @@ subroutine material_init
enddo
endif debugOut
call material_populateGrains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! new mappings
allocate(material_phaseAt(homogenization_maxNgrains,discretization_nElem), source=0)
allocate(material_texture(homogenization_maxNgrains,discretization_nIP,discretization_nElem), source=0) !this is only needed by plasticity nonlocal
allocate(material_EulerAngles(3,homogenization_maxNgrains,discretization_nIP,discretization_nElem),source=0.0_pReal)
do e = 1, discretization_nElem
do i = 1, discretization_nIP
myMicro = discretization_microstructureAt(e)
do c = 1, homogenization_Ngrains(discretization_homogenizationAt(e))
material_phaseAt(c,e) = microstructure_phase(c,myMicro)
material_texture(c,i,e) = microstructure_texture(c,myMicro)
material_EulerAngles(1:3,c,i,e) = texture_Gauss(1:3,material_texture(c,i,e)) ! this is a copy of crystallite_orientation0
enddo
enddo
enddo
deallocate(microstructure_phase)
deallocate(microstructure_texture)
allocate(material_homogenizationAt,source=discretization_homogenizationAt)
allocate(material_homogenizationMemberAt(discretization_nIP,discretization_nElem),source=0)
@ -345,8 +356,6 @@ subroutine material_init
enddo
enddo
allocate(material_phaseAt(homogenization_maxNgrains,discretization_nElem), source=material_phase(:,1,:))
allocate(material_phaseMemberAt(homogenization_maxNgrains,discretization_nIP,discretization_nElem),source=0)
allocate(CounterPhase(size(config_phase)),source=0)
@ -365,8 +374,8 @@ subroutine material_init
#if defined(PETSc) || defined(DAMASK_HDF5)
call results_openJobFile
call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,phase_name)
call results_mapping_materialpoint(material_homogenizationAt,material_homogenizationMemberAt,homogenization_name)
call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,config_name_phase)
call results_mapping_materialpoint(material_homogenizationAt,material_homogenizationMemberAt,config_name_homogenization)
call results_closeJobFile
#endif
@ -375,26 +384,15 @@ subroutine material_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN DEPRECATED
allocate(phaseAt ( homogenization_maxNgrains,discretization_nIP,discretization_nElem),source=0)
allocate(phasememberAt ( homogenization_maxNgrains,discretization_nIP,discretization_nElem),source=0)
allocate(mappingHomogenization (2, discretization_nIP,discretization_nElem),source=0)
allocate(mappingHomogenizationConst( discretization_nIP,discretization_nElem),source=1)
CounterHomogenization=0
CounterPhase =0
do e = 1,discretization_nElem
myHomog = discretization_homogenizationAt(e)
myHomog = discretization_homogenizationAt(e)
do i = 1, discretization_nIP
CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1
mappingHomogenization(1:2,i,e) = [CounterHomogenization(myHomog),huge(1)]
do g = 1,homogenization_Ngrains(myHomog)
myPhase = material_phase(g,i,e)
CounterPhase(myPhase) = CounterPhase(myPhase)+1 ! not distinguishing between instances of same phase
phaseAt(g,i,e) = myPhase
phasememberAt(g,i,e) = CounterPhase(myPhase)
enddo
enddo
enddo
! END DEPRECATED
@ -555,7 +553,7 @@ subroutine material_parseMicrostructure
enddo
enddo
if (dNeq(sum(microstructure_fraction(:,m)),1.0_pReal)) call IO_error(153,ext_msg=microstructure_name(m))
if (dNeq(sum(microstructure_fraction(:,m)),1.0_pReal)) call IO_error(153,ext_msg=config_name_microstructure(m))
enddo
@ -776,41 +774,41 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,&
sizeState,sizeDotState,sizeDeltaState,&
Nslip,Ntwin,Ntrans)
integer, intent(in) :: &
phase, &
NofMyPhase, &
sizeState, &
sizeDotState, &
sizeDeltaState, &
Nslip, &
Ntwin, &
Ntrans
integer, intent(in) :: &
phase, &
NofMyPhase, &
sizeState, &
sizeDotState, &
sizeDeltaState, &
Nslip, &
Ntwin, &
Ntrans
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
plasticState(phase)%Nslip = Nslip
plasticState(phase)%Ntwin = Ntwin
plasticState(phase)%Ntrans= Ntrans
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
plasticState(phase)%Nslip = Nslip
plasticState(phase)%Ntwin = Ntwin
plasticState(phase)%Ntrans= Ntrans
allocate(plasticState(phase)%aTolState (sizeState), source=0.0_pReal)
allocate(plasticState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%aTolState (sizeState), source=0.0_pReal)
allocate(plasticState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 1) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (numerics_integrator == 4) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 5) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 1) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (numerics_integrator == 4) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 5) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
end subroutine material_allocatePlasticState
@ -821,66 +819,35 @@ end subroutine material_allocatePlasticState
subroutine material_allocateSourceState(phase,of,NofMyPhase,&
sizeState,sizeDotState,sizeDeltaState)
integer, intent(in) :: &
phase, &
of, &
NofMyPhase, &
sizeState, sizeDotState,sizeDeltaState
integer, intent(in) :: &
phase, &
of, &
NofMyPhase, &
sizeState, sizeDotState,sizeDeltaState
sourceState(phase)%p(of)%sizeState = sizeState
sourceState(phase)%p(of)%sizeDotState = sizeDotState
sourceState(phase)%p(of)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(of)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
sourceState(phase)%p(of)%sizeState = sizeState
sourceState(phase)%p(of)%sizeDotState = sizeDotState
sourceState(phase)%p(of)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(of)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(sourceState(phase)%p(of)%aTolState (sizeState), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%aTolState (sizeState), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 1) then
allocate(sourceState(phase)%p(of)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (numerics_integrator == 4) &
allocate(sourceState(phase)%p(of)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 5) &
allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 1) then
allocate(sourceState(phase)%p(of)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (numerics_integrator == 4) &
allocate(sourceState(phase)%p(of)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 5) &
allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
end subroutine material_allocateSourceState
!--------------------------------------------------------------------------------------------------
!> @brief populates the grains
!> @details populates the grains by identifying active microstructure/homogenization pairs,
!! calculates the volume of the grains and deals with texture components
!--------------------------------------------------------------------------------------------------
subroutine material_populateGrains
integer :: e,i,c,homog,micro
allocate(material_phase(homogenization_maxNgrains,discretization_nIP,discretization_nElem), source=0)
allocate(material_texture(homogenization_maxNgrains,discretization_nIP,discretization_nElem), source=0)
allocate(material_EulerAngles(3,homogenization_maxNgrains,discretization_nIP,discretization_nElem),source=0.0_pReal)
do e = 1, discretization_nElem
do i = 1, discretization_nIP
homog = discretization_homogenizationAt(e)
micro = discretization_microstructureAt(e)
do c = 1, homogenization_Ngrains(homog)
material_phase(c,i,e) = microstructure_phase(c,micro)
material_texture(c,i,e) = microstructure_texture(c,micro)
material_EulerAngles(1:3,c,i,e) = texture_Gauss(1:3,material_texture(c,i,e))
enddo
enddo
enddo
deallocate(microstructure_phase)
deallocate(microstructure_texture)
end subroutine material_populateGrains
end module material

View File

@ -179,7 +179,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
e = ubound(a,2)
endif
if(present(sortDim)) then
if(present(sortDim)) then
d = sortDim
else
d = 1

View File

@ -291,7 +291,7 @@ program DAMASK_FEM
endif
timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
skipping: if (totalIncsCounter <= restartInc) then ! not yet at restart inc?
skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc?
time = time + timeinc ! just advance time, skip already performed calculation
guess = .true.
else skipping

View File

@ -8,11 +8,10 @@ module mesh
#include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscis.h>
#include <petsc/finclude/petscdmda.h>
use prec
use mesh_base
use PETScdmplex
use PETScdmda
use PETScis
use DAMASK_interface
use IO
use debug
@ -20,6 +19,8 @@ module mesh
use numerics
use FEsolving
use FEM_Zoo
use prec
use mesh_base
implicit none
private
@ -35,13 +36,13 @@ module mesh
mesh_maxNips !< max number of IPs in any CP element
!!!! BEGIN DEPRECATED !!!!!
integer, dimension(:,:), allocatable, public, protected :: &
integer, dimension(:,:), allocatable :: &
mesh_element !DEPRECATED
real(pReal), dimension(:,:), allocatable, public :: &
real(pReal), dimension(:,:), allocatable :: &
mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:), allocatable, public, protected :: &
real(pReal), dimension(:,:), allocatable :: &
mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!)
@ -176,15 +177,13 @@ subroutine mesh_init
endif
enddo
close (FILEUNIT)
endif
if (worldsize > 1) then
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr)
else
call DMClone(globalMesh,geomMesh,ierr)
CHKERRQ(ierr)
else
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr)
endif
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
@ -255,75 +254,66 @@ end function mesh_cellCenterCoordinates
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipVolumes(dimPlex)
PetscInt :: dimPlex
PetscReal :: vol
PetscReal, target :: cent(dimPlex), norm(dimPlex)
PetscReal, pointer :: pCent(:), pNorm(:)
PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr
PetscInt :: dimPlex
PetscReal :: vol
PetscReal, target :: cent(dimPlex), norm(dimPlex)
PetscReal, pointer :: pCent(:), pNorm(:)
PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
endif
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
endif
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
pCent => cent
pNorm => norm
do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr)
mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
enddo
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
pCent => cent
pNorm => norm
do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr)
mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
enddo
end subroutine mesh_FEM_build_ipVolumes
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates'
! Called by all solvers in mesh_init in order to initialize the ip coordinates.
! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus,
! so no need to use this subroutine anymore; Marc however only provides nodal displacements,
! so in this case the ip coordinates are always calculated on the basis of this subroutine.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES,
! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME.
! HAS TO BE CHANGED IN A LATER VERSION.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
qOffset = 0
do qPt = 1, mesh_maxNips
do dirI = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
do dirJ = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0)
enddo
enddo
qOffset = qOffset + dimPlex
enddo
enddo
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
qOffset = 0
do qPt = 1, mesh_maxNips
do dirI = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
do dirJ = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0)
enddo
enddo
qOffset = qOffset + dimPlex
enddo
enddo
end subroutine mesh_FEM_build_ipCoordinates

View File

@ -497,7 +497,7 @@ subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileF
mesh_mapFEtoCPelem(2,cpElem) = cpElem
enddo
call math_sort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2))
call math_sort(mesh_mapFEtoCPelem)
end subroutine mesh_marc_map_elements
@ -532,7 +532,7 @@ subroutine mesh_marc_map_nodes(nNodes,fileUnit)
endif
enddo
620 call math_sort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2))
620 call math_sort(mesh_mapFEtoCPnode)
end subroutine mesh_marc_map_nodes
@ -1262,43 +1262,43 @@ end subroutine mesh_build_ipAreas
!--------------------------------------------------------------------------------------------------
integer function mesh_FEasCP(what,myID)
character(len=*), intent(in) :: what
integer, intent(in) :: myID
character(len=*), intent(in) :: what
integer, intent(in) :: myID
integer, dimension(:,:), pointer :: lookupMap
integer :: lower,upper,center
integer, dimension(:,:), pointer :: lookupMap
integer :: lower,upper,center
mesh_FEasCP = 0
select case(IO_lc(what(1:4)))
case('elem')
lookupMap => mesh_mapFEtoCPelem
case('node')
lookupMap => mesh_mapFEtoCPnode
case default
return
endselect
mesh_FEasCP = 0
select case(IO_lc(what(1:4)))
case('elem')
lookupMap => mesh_mapFEtoCPelem
case('node')
lookupMap => mesh_mapFEtoCPnode
case default
return
endselect
lower = 1
upper = int(size(lookupMap,2),pInt)
lower = 1
upper = int(size(lookupMap,2),pInt)
if (lookupMap(1,lower) == myID) then ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds?
mesh_FEasCP = lookupMap(2,lower)
return
elseif (lookupMap(1,upper) == myID) then
mesh_FEasCP = lookupMap(2,upper)
return
endif
binarySearch: do while (upper-lower > 1)
center = (lower+upper)/2
if (lookupMap(1,center) < myID) then
lower = center
elseif (lookupMap(1,center) > myID) then
upper = center
else
mesh_FEasCP = lookupMap(2,center)
exit
endif
enddo binarySearch
if (lookupMap(1,lower) == myID) then ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds?
mesh_FEasCP = lookupMap(2,lower)
return
elseif (lookupMap(1,upper) == myID) then
mesh_FEasCP = lookupMap(2,upper)
return
endif
binarySearch: do while (upper-lower > 1)
center = (lower+upper)/2
if (lookupMap(1,center) < myID) then
lower = center
elseif (lookupMap(1,center) > myID) then
upper = center
else
mesh_FEasCP = lookupMap(2,center)
exit
endif
enddo binarySearch
end function mesh_FEasCP

View File

@ -13,6 +13,7 @@ module plastic_disloUCLA
use material
use config
use lattice
use discretization
use results
implicit none
@ -295,7 +296,7 @@ subroutine plastic_disloUCLA_init()
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phase == p)
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState

View File

@ -15,9 +15,8 @@ module plastic_dislotwin
use material
use config
use lattice
#if defined(PETSc) || defined(DAMASK_HDF5)
use discretization
use results
#endif
implicit none
private
@ -73,7 +72,7 @@ module plastic_dislotwin
aTol_rho, & !< absolute tolerance for integration of dislocation density
aTol_f_tw, & !< absolute tolerance for integration of twin volume fraction
aTol_f_tr, & !< absolute tolerance for integration of trans volume fraction
gamma_fcc_hex, & !< Free energy difference between austensite and martensite
gamma_fcc_hex, & !< Free energy difference between austensite and martensite
i_tr, & !<
h !< Stack height of hex nucleus
real(pReal), dimension(:), allocatable :: &
@ -89,7 +88,7 @@ module plastic_dislotwin
t_tw, & !< twin thickness [m] for each twin system
CLambdaSlip, & !< Adj. parameter for distance between 2 forest dislocations for each slip system
atomicVolume, &
t_tr, & !< martensite lamellar thickness [m] for each trans system and instance
t_tr, & !< martensite lamellar thickness [m] for each trans system and instance
p, & !< p-exponent in glide velocity
q, & !< q-exponent in glide velocity
r, & !< r-exponent in twin nucleation rate
@ -139,14 +138,14 @@ module plastic_dislotwin
type :: tDislotwinMicrostructure
real(pReal), dimension(:,:), allocatable :: &
Lambda_sl, & !* mean free path between 2 obstacles seen by a moving dislocation
Lambda_tw, & !* mean free path between 2 obstacles seen by a growing twin
Lambda_tr, &!* mean free path between 2 obstacles seen by a growing martensite
Lambda_sl, & !< mean free path between 2 obstacles seen by a moving dislocation
Lambda_tw, & !< mean free path between 2 obstacles seen by a growing twin
Lambda_tr, & !< mean free path between 2 obstacles seen by a growing martensite
tau_pass, &
tau_hat_tw, &
tau_hat_tr, &
f_tw, &
f_tr, &
V_tw, & !< volume of a new twin
V_tr, & !< volume of a new martensite disc
tau_r_tw, & !< stress to bring partials close together (twin)
tau_r_tr !< stress to bring partials close together (trans)
end type tDislotwinMicrostructure
@ -278,7 +277,7 @@ subroutine plastic_dislotwin_init
prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl)
prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl)
prm%v0 = math_expand(prm%v0, prm%N_sl)
prm%b_sl = math_expand(prm%b_sl,prm%N_sl)
prm%b_sl = math_expand(prm%b_sl, prm%N_sl)
prm%Delta_F = math_expand(prm%Delta_F, prm%N_sl)
prm%CLambdaSlip = math_expand(prm%CLambdaSlip, prm%N_sl)
prm%p = math_expand(prm%p, prm%N_sl)
@ -310,23 +309,23 @@ subroutine plastic_dislotwin_init
if (prm%sum_N_tw > 0) then
prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,&
config%getFloats('interaction_twintwin'), &
config%getString('lattice_structure'))
prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,&
config%getFloats('interaction_twintwin'), &
config%getString('lattice_structure'))
prm%b_tw = config%getFloats('twinburgers', requiredSize=size(prm%N_tw))
prm%t_tw = config%getFloats('twinsize', requiredSize=size(prm%N_tw))
prm%r = config%getFloats('r_twin', requiredSize=size(prm%N_tw))
prm%b_tw = config%getFloats('twinburgers', requiredSize=size(prm%N_tw))
prm%t_tw = config%getFloats('twinsize', requiredSize=size(prm%N_tw))
prm%r = config%getFloats('r_twin', requiredSize=size(prm%N_tw))
prm%xc_twin = config%getFloat('xc_twin')
prm%L_tw = config%getFloat('l0_twin')
prm%L_tw = config%getFloat('l0_twin')
prm%i_tw = config%getFloat('cmfptwin')
prm%gamma_char = lattice_characteristicShear_Twin(prm%N_tw,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%C66_tw = lattice_C66_twin(prm%N_tw,prm%C66,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%C66_tw = lattice_C66_twin(prm%N_tw,prm%C66,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if (.not. prm%fccTwinTransNucleation) then
prm%dot_N_0_tw = config%getFloats('ndot0_twin')
@ -339,14 +338,16 @@ subroutine plastic_dislotwin_init
prm%r = math_expand(prm%r,prm%N_tw)
else
allocate(prm%t_tw(0))
allocate(prm%b_tw(0))
allocate(prm%r(0))
allocate(prm%gamma_char(0))
allocate(prm%t_tw (0))
allocate(prm%b_tw (0))
allocate(prm%r (0))
allocate(prm%h_tw_tw (0,0))
endif
!--------------------------------------------------------------------------------------------------
! transformation related parameters
prm%N_tr = config%getInts('ntrans', defaultVal=emptyIntArray)
prm%N_tr = config%getInts('ntrans', defaultVal=emptyIntArray)
prm%sum_N_tr = sum(prm%N_tr)
if (prm%sum_N_tr > 0) then
prm%b_tr = config%getFloats('transburgers')
@ -383,8 +384,10 @@ subroutine plastic_dislotwin_init
prm%s = config%getFloats('s_trans',defaultVal=[0.0_pReal])
prm%s = math_expand(prm%s,prm%N_tr)
else
allocate(prm%t_tr(0))
allocate(prm%b_tr(0))
allocate(prm%t_tr (0))
allocate(prm%b_tr (0))
allocate(prm%s (0))
allocate(prm%h_tr_tr(0,0))
endif
if (sum(prm%N_tw) > 0 .or. prm%sum_N_tr > 0) then
@ -452,42 +455,33 @@ subroutine plastic_dislotwin_init
do i= 1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('edge_density')
case ('rho_mob')
outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl > 0)
outputSize = prm%sum_N_sl
case ('dipole_density')
case ('rho_dip')
outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl > 0)
outputSize = prm%sum_N_sl
case ('shear_rate_slip','shearrate_slip')
outputID = merge(dot_gamma_sl_ID,undefined_ID,prm%sum_N_sl > 0)
outputSize = prm%sum_N_sl
case ('accumulated_shear_slip')
case ('gamma_sl')
outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl > 0)
outputSize = prm%sum_N_sl
case ('mfp_slip')
case ('lambda_sl')
outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl > 0)
outputSize = prm%sum_N_sl
case ('resolved_stress_slip')
outputID = merge(resolved_stress_slip_ID,undefined_ID,prm%sum_N_sl > 0)
outputSize = prm%sum_N_sl
case ('threshold_stress_slip')
case ('tau_pass')
outputID= merge(threshold_stress_slip_ID,undefined_ID,prm%sum_N_sl > 0)
outputSize = prm%sum_N_sl
case ('twin_fraction')
case ('f_tw')
outputID = merge(f_tw_ID,undefined_ID,prm%sum_N_tw >0)
outputSize = prm%sum_N_tw
case ('mfp_twin')
case ('lambda_tw')
outputID = merge(Lambda_tw_ID,undefined_ID,prm%sum_N_tw >0)
outputSize = prm%sum_N_tw
case ('resolved_stress_twin')
outputID = merge(resolved_stress_twin_ID,undefined_ID,prm%sum_N_tw >0)
outputSize = prm%sum_N_tw
case ('threshold_stress_twin')
case ('tau_hat_tw')
outputID = merge(tau_hat_tw_ID,undefined_ID,prm%sum_N_tw >0)
outputSize = prm%sum_N_tw
case ('strain_trans_fraction')
case ('f_tr')
outputID = f_tr_ID
outputSize = prm%sum_N_tr
@ -503,7 +497,7 @@ subroutine plastic_dislotwin_init
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phase == p)
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
+ size(['f_tw']) * prm%sum_N_tw &
+ size(['f_tr']) * prm%sum_N_tr
@ -551,18 +545,18 @@ subroutine plastic_dislotwin_init
dot%f_tr=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tr
allocate(dst%Lambda_sl (prm%sum_N_sl, NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_pass (prm%sum_N_sl, NipcMyPhase),source=0.0_pReal)
allocate(dst%Lambda_sl (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_pass (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal)
allocate(dst%Lambda_tw (prm%sum_N_tw, NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_hat_tw (prm%sum_N_tw, NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_r_tw (prm%sum_N_tw, NipcMyPhase),source=0.0_pReal)
allocate(dst%f_tw (prm%sum_N_tw, NipcMyPhase),source=0.0_pReal)
allocate(dst%Lambda_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_hat_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_r_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal)
allocate(dst%V_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal)
allocate(dst%Lambda_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_hat_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal)
allocate(dst%tau_r_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal)
allocate(dst%f_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal)
allocate(dst%V_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
@ -590,9 +584,9 @@ function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
of
real(pReal) :: f_unrotated
of = phasememberAt(ipc,ip,el)
associate(prm => param(phase_plasticityInstance(material_phase(ipc,ip,el))),&
stt => state(phase_plasticityInstance(material_phase(ipc,ip,el))))
of = material_phasememberAt(ipc,ip,el)
associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),&
stt => state(phase_plasticityInstance(material_phaseAT(ipc,el))))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,of)) &
@ -742,10 +736,12 @@ subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
of
integer :: i
real(pReal) :: f_unrotated,&
VacancyDiffusion,&
rho_dip_distance, ClimbVelocity, &
tau
real(pReal) :: &
f_unrotated, &
VacancyDiffusion, &
rho_dip_distance, &
v_cl, &
tau
real(pReal), dimension(param(instance)%sum_N_sl) :: &
dot_rho_dip_formation, &
dot_rho_dip_climb, &
@ -757,7 +753,7 @@ subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
dot_gamma_tr
associate(prm => param(instance), stt => state(instance), &
dot => dotstate(instance), dst => dependentState(instance))
dot => dotState(instance), dst => dependentState(instance))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,of)) &
@ -790,9 +786,9 @@ subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
if (dEq0(rho_dip_distance-rho_dip_distance_min(i))) then
dot_rho_dip_climb(i) = 0.0_pReal
else
ClimbVelocity = 3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume(i) &
/ (2.0_pReal*PI*kB*T*(rho_dip_distance+rho_dip_distance_min(i)))
dot_rho_dip_climb(i) = 4.0_pReal*ClimbVelocity*stt%rho_dip(i,of) &
v_cl = 3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume(i) &
/ (2.0_pReal*PI*kB*T*(rho_dip_distance+rho_dip_distance_min(i)))
dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) &
/ (rho_dip_distance-rho_dip_distance_min(i))
endif
endif significantSlipStress
@ -811,7 +807,7 @@ subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
dot%f_tw(:,of) = f_unrotated*dot_gamma_twin/prm%gamma_char
call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr)
dot%f_tw(:,of) = f_unrotated*dot_gamma_tr
dot%f_tr(:,of) = f_unrotated*dot_gamma_tr
end associate
@ -834,18 +830,17 @@ subroutine plastic_dislotwin_dependentState(T,instance,of)
real(pReal) :: &
sumf_twin,SFE,sumf_trans
real(pReal), dimension(param(instance)%sum_N_sl) :: &
inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation
inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation
inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation
inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation
real(pReal), dimension(param(instance)%sum_N_tw) :: &
inv_lambda_tw_tw !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
f_over_t_tw
real(pReal), dimension(param(instance)%sum_N_tr) :: &
inv_lambda_tr_tr !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite
real(pReal), dimension(:), allocatable :: &
x0, &
f_over_t_tw, &
inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite
f_over_t_tr
real(pReal), dimension(:), allocatable :: &
x0
associate(prm => param(instance),&
@ -858,9 +853,9 @@ subroutine plastic_dislotwin_dependentState(T,instance,of)
SFE = prm%SFE_0K + prm%dSFE_dT * T
!* rescaled volume fraction for topology
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw !ToDo: this is per system
f_over_t_tr = sumf_trans/prm%t_tr !ToDo: But this not ...
!Todo: Physically ok, but naming could be adjusted
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ...
f_over_t_tr = sumf_trans/prm%t_tr ! but this not
! ToDo ...Physically correct, but naming could be adjusted
forall (i = 1:prm%sum_N_sl) &
@ -872,30 +867,22 @@ subroutine plastic_dislotwin_dependentState(T,instance,of)
if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
!ToDo: needed? if (prm%sum_N_tw > 0) &
inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans)
!ToDo: needed? if (prm%sum_N_tr > 0) &
inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans)
if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: Change order
dst%Lambda_sl(:,of) = &
prm%D/(1.0_pReal+prm%D*&
(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr))
else
dst%Lambda_sl(:,of) = prm%D &
/ (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct?
endif
if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: better logic needed here
dst%Lambda_sl(:,of) = prm%D &
/ (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr))
else
dst%Lambda_sl(:,of) = prm%D &
/ (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct?
endif
dst%Lambda_tw(:,of) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
@ -906,16 +893,16 @@ subroutine plastic_dislotwin_dependentState(T,instance,of)
!* threshold stress for growing twin/martensite
if(prm%sum_N_tw == prm%sum_N_sl) &
dst%tau_hat_tw(:,of) = &
(SFE/(3.0_pReal*prm%b_tw)+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl)) ! slip burgers here correct?
dst%tau_hat_tw(:,of) = SFE/(3.0_pReal*prm%b_tw) &
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip burgers here correct?
if(prm%sum_N_tr == prm%sum_N_sl) &
dst%tau_hat_tr(:,of) = &
(SFE/(3.0_pReal*prm%b_tr) + 3.0_pReal*prm%b_tr*prm%mu/&
(prm%L_tr*prm%b_sl) + prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) )
dst%tau_hat_tr(:,of) = SFE/(3.0_pReal*prm%b_tr) &
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip burgers here correct?
+ prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr)
dst%f_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal
dst%f_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal
dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal
dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal
x0 = prm%mu*prm%b_tw**2.0_pReal/(SFE*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip and is the same for twin and trans
@ -1174,12 +1161,11 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,&
isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < dst%tau_r_tw(i,of)) then
if (tau(i) < dst%tau_r_tw(i,of)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state
(prm%L_tw*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*&
(dst%tau_r_tw(i,of)-tau)))
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,of)-tau(i)))) ! P_ncs
else
Ndot0=0.0_pReal
end if
@ -1189,8 +1175,8 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,&
enddo
significantStress: where(tau > tol_math_check)
StressRatio_r = (dst%tau_hat_tw(:,of)/tau)**prm%r
dot_gamma_twin = prm%gamma_char * dst%f_tw(:,of) * Ndot0*exp(-StressRatio_r)
StressRatio_r = (dst%tau_hat_tw(:,of)/tau)**prm%r
dot_gamma_twin = prm%gamma_char * dst%V_tw(:,of) * Ndot0*exp(-StressRatio_r)
ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r
else where significantStress
dot_gamma_twin = 0.0_pReal
@ -1232,7 +1218,6 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
ddot_gamma_dtau
integer :: i,s1,s2
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
do i = 1, prm%sum_N_tr
@ -1240,12 +1225,11 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < dst%tau_r_tr(i,of)) then
if (tau(i) < dst%tau_r_tr(i,of)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state
(prm%L_tr*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*&
(dst%tau_r_tr(i,of)-tau)))
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,of)-tau(i)))) ! P_ncs
else
Ndot0=0.0_pReal
end if
@ -1255,9 +1239,9 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
enddo
significantStress: where(tau > tol_math_check)
StressRatio_s = (dst%tau_hat_tr(:,of)/tau)**prm%s
dot_gamma_tr = dst%f_tr(:,of) * Ndot0*exp(-StressRatio_s)
ddot_gamma_dtau = (dot_gamma_tr*prm%r/tau)*StressRatio_s
StressRatio_s = (dst%tau_hat_tr(:,of)/tau)**prm%s
dot_gamma_tr = dst%V_tr(:,of) * Ndot0*exp(-StressRatio_s)
ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s
else where significantStress
dot_gamma_tr = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal

View File

@ -14,9 +14,8 @@ module plastic_isotropic
use IO
use material
use config
#if defined(PETSc) || defined(DAMASK_HDF5)
use discretization
use results
#endif
implicit none
private
@ -127,8 +126,8 @@ subroutine plastic_isotropic_init
config => config_phase(p))
#ifdef DEBUG
if (p==material_phase(debug_g,debug_i,debug_e)) &
prm%of_debug = phasememberAt(debug_g,debug_i,debug_e)
if (p==material_phaseAt(debug_g,debug_e)) &
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
#endif
prm%xi_0 = config%getFloat('tau0')
@ -190,7 +189,7 @@ subroutine plastic_isotropic_init
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phase == p)
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['xi ','accumulated_shear'])
sizeState = sizeDotState

View File

@ -13,9 +13,8 @@ module plastic_kinehardening
use material
use config
use lattice
#if defined(PETSc) || defined(DAMASK_HDF5)
use discretization
use results
#endif
implicit none
private
@ -146,8 +145,8 @@ subroutine plastic_kinehardening_init
config => config_phase(p))
#ifdef DEBUG
if (p==material_phase(debug_g,debug_i,debug_e)) then
prm%of_debug = phasememberAt(debug_g,debug_i,debug_e)
if (p==material_phaseAt(debug_g,debug_e)) then
prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
endif
#endif
@ -257,7 +256,7 @@ subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phase == p)
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip
sizeState = sizeDotState + sizeDeltaState

View File

@ -6,6 +6,7 @@
!--------------------------------------------------------------------------------------------------
module plastic_none
use material
use discretization
use debug
implicit none
@ -36,7 +37,7 @@ subroutine plastic_none_init
do p = 1, size(phase_plasticity)
if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
NipcMyPhase = count(material_phase == p)
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
call material_allocatePlasticState(p,NipcMyPhase,0,0,0, &
0,0,0)
plasticState(p)%sizePostResults = 0

View File

@ -556,7 +556,7 @@ subroutine plastic_nonlocal_init
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NofMyPhase=count(material_phase==p)
NofMyPhase = count(material_phaseAt==p) * discretization_nIP
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@ -677,7 +677,7 @@ subroutine plastic_nonlocal_init
allocate(iD(maxval(totalNslip),2,maxNinstances), source=0)
initializeInstances: do p = 1, size(phase_plasticity)
NofMyPhase=count(material_phase==p)
NofMyPhase = count(material_phaseAt==p) * discretization_nIP
myPhase2: if (phase_plasticity(p) == PLASTICITY_NONLOCAL_ID) then
!*** determine indices to state array
@ -766,7 +766,7 @@ subroutine plastic_nonlocal_init
! get the total volume of the instance
do e = 1,discretization_nElem
do i = 1,discretization_nIP
if (material_phase(1,i,e) == phase) volume(phasememberAt(1,i,e)) = IPvolume(i,e)
if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e)
enddo
enddo
totalVolume = sum(volume)
@ -854,29 +854,29 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
invConnections
real(pReal), dimension(3,nIPneighbors) :: &
connection_latticeConf
real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
rhoExcess
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
rho_edg_delta, &
rho_scr_delta
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
rho, &
rho_neighbor
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))), &
totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))), &
totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
myInteractionMatrix ! corrected slip interaction matrix
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),nIPneighbors) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),nIPneighbors) :: &
rho_edg_delta_neighbor, &
rho_scr_delta_neighbor
real(pReal), dimension(2,maxval(totalNslip),nIPneighbors) :: &
neighbor_rhoExcess, & ! excess density at neighboring material point
neighbor_rhoTotal ! total density at neighboring material point
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),2) :: &
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: &
m ! direction of dislocation motion
ph = phaseAt(1,ip,el)
of = phasememberAt(1,ip,el)
ph = material_phaseAt(1,el)
of = material_phasememberAt(1,ip,el)
instance = phase_plasticityInstance(ph)
associate(prm => param(instance),dst => microstructure(instance), stt => state(instance))
@ -935,9 +935,9 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
no = phasememberAt(1,neighbor_ip,neighbor_el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
if (neighbor_el > 0 .and. neighbor_ip > 0) then
neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el))
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
if (neighbor_instance == instance) then
nRealNeighbors = nRealNeighbors + 1.0_pReal
@ -1202,22 +1202,22 @@ subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, &
of, & !offset
t, & !< dislocation type
s !< index of my current slip system
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),8) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: &
rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
rho
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
v, & !< velocity
tauNS, & !< resolved shear stress including non Schmid and backstress terms
dv_dtau, & !< velocity derivative with respect to the shear stress
dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate
!*** shortcut for mapping
ph = phaseAt(1,ip,el)
of = phasememberAt(1,ip,el)
ph = material_phaseAt(1,el)
of = material_phasememberAt(1,ip,el)
instance = phase_plasticityInstance(ph)
associate(prm => param(instance),dst=>microstructure(instance))
@ -1323,23 +1323,23 @@ subroutine plastic_nonlocal_deltaState(Mp,ip,el)
c, & ! character of dislocation
t, & ! type of dislocation
s ! index of my current slip system
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
deltaRhoRemobilization, & ! density increment by remobilization
deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
rho ! current dislocation densities
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
v ! dislocation glide velocity
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
tau ! current resolved shear stress
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),2) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: &
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
dUpper, & ! current maximum stable dipole distance for edges and screws
dUpperOld, & ! old maximum stable dipole distance for edges and screws
deltaDUpper ! change in maximum stable dipole distance for edges and screws
ph = phaseAt(1,ip,el)
of = phasememberAt(1,ip,el)
ph = material_phaseAt(1,el)
of = material_phasememberAt(1,ip,el)
instance = phase_plasticityInstance(ph)
associate(prm => param(instance),dst => microstructure(instance),del => deltaState(instance))
@ -1459,7 +1459,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
np,& !< neighbour phase shortcut
topp, & !< type of dislocation with opposite sign to t
s !< index of my current slip system
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
rho, &
rhoDot, & !< density evolution
rhoDotMultiplication, & !< density evolution by multiplication
@ -1467,24 +1467,24 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),8) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
v, & !< current dislocation glide velocity
my_v, & !< dislocation glide velocity of central ip
neighbor_v, & !< dislocation glide velocity of enighboring ip
gdot !< shear rates
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
tau, & !< current resolved shear stress
vClimb !< climb velocity of edge dipoles
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),2) :: &
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: &
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: &
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
m !< direction of dislocation motion
real(pReal), dimension(3,3) :: &
my_F, & !< my total deformation gradient
@ -1507,15 +1507,15 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
considerEnteringFlux, &
considerLeavingFlux
p = phaseAt(1,ip,el)
o = phasememberAt(1,ip,el)
p = material_phaseAt(1,el)
o = material_phasememberAt(1,ip,el)
if (timestep <= 0.0_pReal) then
plasticState(p)%dotState = 0.0_pReal
return
endif
ph = material_phase(1,ip,el)
ph = material_phaseAt(1,el)
instance = phase_plasticityInstance(ph)
associate(prm => param(instance),dst => microstructure(instance),dot => dotState(instance),stt => state(instance))
ns = totalNslip(instance)
@ -1592,7 +1592,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
!****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal plasticity)
rhoDotFlux = 0.0_pReal
if (.not. phase_localPlasticity(material_phase(1,ip,el))) then
if (.not. phase_localPlasticity(material_phaseAt(1,el))) then
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(gdot) > 0.0_pReal & ! any active slip system ...
@ -1630,8 +1630,8 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
neighbor_n = IPneighborhood(3,n,ip,el)
np = phaseAt(1,neighbor_ip,neighbor_el)
no = phasememberAt(1,neighbor_ip,neighbor_el)
np = material_phaseAt(1,neighbor_el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
opposite_neighbor = n + mod(n,2) - mod(n+1,2)
opposite_el = IPneighborhood(1,opposite_neighbor,ip,el)
@ -1639,7 +1639,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el))
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el)
neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el))
Favg = 0.5_pReal * (my_F + neighbor_F)
@ -1661,7 +1661,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below
neighbor_rhoSgl = 0.0_pReal
if (neighbor_n > 0) then
if (phase_plasticity(material_phase(1,neighbor_ip,neighbor_el)) == PLASTICITY_NONLOCAL_ID &
if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID &
.and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) &
considerEnteringFlux = .true.
endif
@ -1714,7 +1714,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
considerLeavingFlux = .true.
if (opposite_n > 0) then
if (phase_plasticity(material_phase(1,opposite_ip,opposite_el)) /= PLASTICITY_NONLOCAL_ID) &
if (phase_plasticity(material_phaseAt(1,opposite_el)) /= PLASTICITY_NONLOCAL_ID) &
considerLeavingFlux = .false.
endif
@ -1905,20 +1905,20 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
s2 ! slip system index (my neighbor)
real(pReal), dimension(4) :: &
absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor
real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),&
totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),&
nIPneighbors) :: &
my_compatibility ! my_compatibility for current element and ip
real(pReal) :: &
my_compatibilitySum, &
thresholdValue, &
nThresholdValues
logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: &
logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: &
belowThreshold
type(rotation) :: rot
Nneighbors = nIPneighbors
ph = material_phase(1,i,e)
ph = material_phaseAt(1,e)
textureID = material_texture(1,i,e)
instance = phase_plasticityInstance(ph)
ns = totalNslip(instance)
@ -1950,7 +1950,7 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
!* we consider this to be a real "physical" phase boundary, so completely incompatible.
!* If one of the two phases has a local plasticity law,
!* we do not consider this to be a phase boundary, so completely compatible.
neighbor_phase = material_phase(1,neighbor_i,neighbor_e)
neighbor_phase = material_phaseAt(1,neighbor_e)
if (neighbor_phase /= ph) then
if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph))&
forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,n) = 0.0_pReal

View File

@ -12,9 +12,8 @@ module plastic_phenopowerlaw
use material
use config
use lattice
#if defined(PETSc) || defined(DAMASK_HDF5)
use discretization
use results
#endif
implicit none
private
@ -314,7 +313,7 @@ subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phase == p)
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip &
+ size(['tau_twin ','gamma_twin']) * prm%totalNtwin
sizeState = sizeDotState

View File

@ -7,43 +7,42 @@
!--------------------------------------------------------------------------------------------------
subroutine quit(stop_id)
#include <petsc/finclude/petscsys.h>
use PetscSys
#ifdef _OPENMP
use MPI, only: &
MPI_finalize
use MPI
#endif
use PetscSys
use hdf5
use hdf5
implicit none
integer, intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
integer :: error
PetscErrorCode :: ierr = 0
implicit none
integer, intent(in) :: stop_id
integer, dimension(8) :: dateAndTime
integer :: error
PetscErrorCode :: ierr = 0
call h5open_f(error)
if (error /= 0) write(6,'(a,i5)') ' Error in h5open_f ',error ! prevents error if not opened yet
call h5close_f(error)
if (error /= 0) write(6,'(a,i5)') ' Error in h5close_f ',error
call h5open_f(error)
if (error /= 0) write(6,'(a,i5)') ' Error in h5open_f ',error ! prevents error if not opened yet
call h5close_f(error)
if (error /= 0) write(6,'(a,i5)') ' Error in h5close_f ',error
call PETScFinalize(ierr)
CHKERRQ(ierr)
call PETScFinalize(ierr)
CHKERRQ(ierr)
#ifdef _OPENMP
call MPI_finalize(error)
if (error /= 0) write(6,'(a,i5)') ' Error in MPI_finalize',error
call MPI_finalize(error)
if (error /= 0) write(6,'(a,i5)') ' Error in MPI_finalize',error
#endif
call date_and_time(values = dateAndTime)
write(6,'(/,a)') ' DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
call date_and_time(values = dateAndTime)
write(6,'(/,a)') ' DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
if (stop_id == 0 .and. ierr == 0 .and. error == 0) stop 0 ! normal termination
if (stop_id == 2 .and. ierr == 0 .and. error == 0) stop 2 ! not all incs converged
stop 1 ! error (message from IO_error)
if (stop_id == 0 .and. ierr == 0 .and. error == 0) stop 0 ! normal termination
if (stop_id == 2 .and. ierr == 0 .and. error == 0) stop 2 ! not all incs converged
stop 1 ! error (message from IO_error)
end subroutine quit

View File

@ -17,8 +17,7 @@ module results
private
#if defined(PETSc) || defined(DAMASK_HDF5)
integer(HID_T), public, protected :: tempCoordinates, tempResults
integer(HID_T), private :: resultsFile, currentIncID, plist_id
integer(HID_T) :: resultsFile
interface results_writeDataset

View File

@ -10,6 +10,7 @@ module source_damage_anisoBrittle
use IO
use math
use material
use discretization
use config
use lattice
@ -164,7 +165,7 @@ subroutine source_damage_anisoBrittle_init
end associate
phase = p
NofMyPhase=count(material_phase==phase)
NofMyPhase=count(material_phaseAt==phase) * discretization_nIP
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
@ -202,8 +203,8 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
phase = material_phaseAt(ipc,el)
constituent = material_phasememberAt(ipc,ip,el)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
homog = material_homogenizationAt(el)

View File

@ -5,55 +5,56 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_anisoDuctile
use prec
use debug
use IO
use math
use material
use config
use prec
use debug
use IO
use math
use discretization
use material
use config
implicit none
private
implicit none
private
integer, dimension(:), allocatable, public, protected :: &
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism
integer, dimension(:), allocatable, public, protected :: &
source_damage_anisoDuctile_offset, & !< which source is my current damage 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
integer, dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_output !< name of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_output !< name of each post result output
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
aTol, &
N
real(pReal), dimension(:), allocatable :: &
critPlasticStrain
integer :: &
totalNslip
integer, dimension(:), allocatable :: &
Nslip
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
end type tParameters
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
aTol, &
N
real(pReal), dimension(:), allocatable :: &
critPlasticStrain
integer :: &
totalNslip
integer, dimension(:), allocatable :: &
Nslip
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_damage_anisoDuctile_init, &
source_damage_anisoDuctile_dotState, &
source_damage_anisoDuctile_getRateAndItsTangent, &
source_damage_anisoDuctile_postResults
public :: &
source_damage_anisoDuctile_init, &
source_damage_anisoDuctile_dotState, &
source_damage_anisoDuctile_getRateAndItsTangent, &
source_damage_anisoDuctile_postResults
contains
@ -64,193 +65,196 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_init
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p ,i
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p ,i
integer, dimension(0), parameter :: emptyIntArray = [integer::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
outputID
integer, dimension(0), parameter :: emptyIntArray = [integer::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
character(len=pStringLen) :: &
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'
Ninstance = count(phase_source == SOURCE_damage_anisoDuctile_ID)
if (Ninstance == 0) return
Ninstance = count(phase_source == SOURCE_damage_anisoDuctile_ID)
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoDuctile_offset(size(config_phase)), source=0)
allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0)
do phase = 1, size(config_phase)
source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_damage_anisoDuctile_ID) &
source_damage_anisoDuctile_offset(phase) = source
enddo
enddo
allocate(source_damage_anisoDuctile_offset(size(config_phase)), source=0)
allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0)
do phase = 1, size(config_phase)
source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_damage_anisoDuctile_ID) &
source_damage_anisoDuctile_offset(phase) = source
enddo
enddo
allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoDuctile_output = ''
allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoDuctile_output = ''
allocate(param(Ninstance))
allocate(param(Ninstance))
do p=1, size(config_phase)
if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISODUCTILE_ID)) cycle
associate(prm => param(source_damage_anisoDuctile_instance(p)), &
config => config_phase(p))
do p=1, size(config_phase)
if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISODUCTILE_ID)) cycle
associate(prm => param(source_damage_anisoDuctile_instance(p)), &
config => config_phase(p))
prm%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal)
prm%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal)
prm%N = config%getFloat('anisoductile_ratesensitivity')
prm%totalNslip = sum(prm%Nslip)
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol'
prm%N = config%getFloat('anisoductile_ratesensitivity')
prm%totalNslip = sum(prm%Nslip)
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity'
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip))
prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip))
! expand: family => system
prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip)
! expand: family => system
prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip)
if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain'
if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')')
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
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)
prm%outputID = [prm%outputID, damage_drivingforce_ID]
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)
prm%outputID = [prm%outputID, damage_drivingforce_ID]
end select
end select
enddo
enddo
end associate
end associate
phase = p
phase = p
NofMyPhase=count(material_phase==phase)
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
NofMyPhase=count(material_phaseAt==phase) * discretization_nIP
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
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
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
enddo
enddo
end subroutine source_damage_anisoDuctile_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: &
phase, &
constituent, &
sourceOffset, &
homog, damageOffset, &
instance, &
f, i
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: &
phase, &
constituent, &
sourceOffset, &
homog, damageOffset, &
instance, &
i
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
phase = material_phaseAt(ipc,el)
constituent = material_phasememberAt(ipc,ip,el)
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
do i = 1, param(instance)%totalNslip
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
plasticState(phase)%slipRate(i,constituent)/ &
((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i)
enddo
do i = 1, param(instance)%totalNslip
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
plasticState(phase)%slipRate(i,constituent)/ &
((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i)
enddo
end subroutine source_damage_anisoDuctile_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, &
constituent
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
sourceOffset
integer, intent(in) :: &
phase, &
constituent
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
sourceOffset
sourceOffset = source_damage_anisoDuctile_offset(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
localphiDot = 1.0_pReal &
- sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi
localphiDot = 1.0_pReal &
- sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
end subroutine source_damage_anisoDuctile_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results
!--------------------------------------------------------------------------------------------------
function source_damage_anisoDuctile_postResults(phase, constituent)
integer, intent(in) :: &
phase, &
constituent
real(pReal), dimension(sum(source_damage_anisoDuctile_sizePostResult(:, &
source_damage_anisoDuctile_instance(phase)))) :: &
source_damage_anisoDuctile_postResults
integer, intent(in) :: &
phase, &
constituent
real(pReal), dimension(sum(source_damage_anisoDuctile_sizePostResult(:, &
source_damage_anisoDuctile_instance(phase)))) :: &
source_damage_anisoDuctile_postResults
integer :: &
instance, sourceOffset, o, c
integer :: &
instance, sourceOffset, o, c
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
c = 0
c = 0
do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID)
source_damage_anisoDuctile_postResults(c+1) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID)
source_damage_anisoDuctile_postResults(c+1) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
end select
enddo
end select
enddo
end function source_damage_anisoDuctile_postResults
end module source_damage_anisoDuctile

View File

@ -9,6 +9,7 @@ module source_damage_isoBrittle
use debug
use IO
use math
use discretization
use material
use config
@ -133,7 +134,7 @@ subroutine source_damage_isoBrittle_init
phase = p
NofMyPhase=count(material_phase==phase)
NofMyPhase = count(material_phaseAt==phase) * discretization_nIP
instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase)
@ -164,8 +165,8 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
strain(6), &
strainenergy
phase = phaseAt(ipc,ip,el) !< phase ID at ipc,ip,el
constituent = phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el
phase = material_phaseAt(ipc,el) !< phase ID at ipc,ip,el
constituent = material_phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el
! ToDo: capability for multiple instances of SAME source within given phase. Needs Ninstance loop from here on!
instance = source_damage_isoBrittle_instance(phase) !< instance of damage_isoBrittle source
sourceOffset = source_damage_isoBrittle_offset(phase)

View File

@ -8,6 +8,7 @@ module source_damage_isoDuctile
use prec
use debug
use IO
use discretization
use material
use config
@ -132,7 +133,7 @@ subroutine source_damage_isoDuctile_init
end associate
phase = p
NofMyPhase=count(material_phase==phase)
NofMyPhase=count(material_phaseAt==phase) * discretization_nIP
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
@ -157,8 +158,8 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
integer :: &
phase, constituent, instance, homog, sourceOffset, damageOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
phase = material_phaseAt(ipc,el)
constituent = material_phasememberAt(ipc,ip,el)
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
homog = material_homogenizationAt(el)

Some files were not shown because too many files have changed in this diff Show More