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

Merge remote-tracking branch 'origin/development' into no-crystallite-dPdF

This commit is contained in:
Martin Diehl 2020-10-03 08:34:02 +02:00
parent 57174d0aba c9ead8bf59
commit cdf3c8cdee
31 changed files with 387473 additions and 509 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 25ce39dd0f5bf49dc5c2bec20767d93d2b76d353 Subproject commit 3b498f5cb3c50e669588106de1b4cdc4c03ffff1

2
VERSION
View File

@ -1 +1 @@
v3.0.0-alpha-374-gc545ad869 v3.0.0-alpha-433-g58229b885

207
examples/SpectralMethod/Polycrystal/material.yaml
View File

@ -1,107 +1,109 @@
homogenization: homogenization:
SX: SX:
mech: {type: none} mech: {type: none}
microstructure:
- constituents: material:
- fraction: 1.0 - constituents:
orientation: [1.0, 0.0, 0.0, 0.0] - fraction: 1.0
phase: Aluminum O: [1.0, 0.0, 0.0, 0.0]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434] - fraction: 1.0
phase: Aluminum O: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945] - fraction: 1.0
phase: Aluminum O: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106] - fraction: 1.0
phase: Aluminum O: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.33012772942625784, -0.6781865350268957, 0.6494525351030648, 0.09638521992649676] - fraction: 1.0
phase: Aluminum O: [0.33012772942625784, -0.6781865350268957, 0.6494525351030648, 0.09638521992649676]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.43596817439583935, -0.5982537129781701, 0.046599032277502436, 0.6707106499919265] - fraction: 1.0
phase: Aluminum O: [0.43596817439583935, -0.5982537129781701, 0.046599032277502436, 0.6707106499919265]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.169734823419553, -0.699615227367322, -0.6059581215838098, -0.33844257746495854] - fraction: 1.0
phase: Aluminum O: [0.169734823419553, -0.699615227367322, -0.6059581215838098, -0.33844257746495854]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.9698864809294915, 0.1729052643205874, -0.15948307917616958, 0.06315956884687175] - fraction: 1.0
phase: Aluminum O: [0.9698864809294915, 0.1729052643205874, -0.15948307917616958, 0.06315956884687175]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.46205660912967883, 0.3105054068891252, -0.617849551030653, 0.555294529545738] - fraction: 1.0
phase: Aluminum O: [0.46205660912967883, 0.3105054068891252, -0.617849551030653, 0.555294529545738]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.4512443497461787, -0.7636045534540555, -0.04739348426715133, -0.45939142396805815] - fraction: 1.0
phase: Aluminum O: [0.4512443497461787, -0.7636045534540555, -0.04739348426715133, -0.45939142396805815]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.2161856212656443, -0.6581450184826598, -0.5498086209601588, 0.4667112513346289] - fraction: 1.0
phase: Aluminum O: [0.2161856212656443, -0.6581450184826598, -0.5498086209601588, 0.4667112513346289]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.8753220715350803, -0.4561599367657419, -0.13298279533852678, -0.08969369719975541] - fraction: 1.0
phase: Aluminum O: [0.8753220715350803, -0.4561599367657419, -0.13298279533852678, -0.08969369719975541]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.11908260752431069, 0.18266024809834172, -0.7144822594012615, -0.664807992845101] - fraction: 1.0
phase: Aluminum O: [0.11908260752431069, 0.18266024809834172, -0.7144822594012615, -0.664807992845101]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.751104669484278, 0.5585633382623958, -0.34579336397009175, 0.06538900566860861] - fraction: 1.0
phase: Aluminum O: [0.751104669484278, 0.5585633382623958, -0.34579336397009175, 0.06538900566860861]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.08740438971703973, 0.8991264096610437, -0.4156704205935976, 0.10559485570696363] - fraction: 1.0
phase: Aluminum O: [0.08740438971703973, 0.8991264096610437, -0.4156704205935976, 0.10559485570696363]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.5584325870096193, 0.6016408353068798, -0.14280340445801173, 0.5529814994483859] - fraction: 1.0
phase: Aluminum O: [0.5584325870096193, 0.6016408353068798, -0.14280340445801173, 0.5529814994483859]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.4052725440888093, 0.25253073423599154, 0.5693263597910454, -0.669215876471182] - fraction: 1.0
phase: Aluminum O: [0.4052725440888093, 0.25253073423599154, 0.5693263597910454, -0.669215876471182]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.7570164606888676, 0.15265448024694664, -0.5998021466848317, 0.20942796551297105] - fraction: 1.0
phase: Aluminum O: [0.7570164606888676, 0.15265448024694664, -0.5998021466848317, 0.20942796551297105]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.6987659297138081, -0.132172211261028, -0.19693254724422338, 0.6748883269678543] - fraction: 1.0
phase: Aluminum O: [0.6987659297138081, -0.132172211261028, -0.19693254724422338, 0.6748883269678543]
homogenization: SX phase: Aluminum
- constituents: homogenization: SX
- fraction: 1.0 - constituents:
orientation: [0.7729330445886478, 0.21682179052722322, -0.5207379472917645, 0.2905078484066341] - fraction: 1.0
phase: Aluminum O: [0.7729330445886478, 0.21682179052722322, -0.5207379472917645, 0.2905078484066341]
homogenization: SX phase: Aluminum
homogenization: SX
phase: phase:
Aluminum: Aluminum:
elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke} elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}
@ -117,7 +119,6 @@ phase:
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4] h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
n_sl: 20 n_sl: 20
output: [xi_sl] output: [xi_sl]
type: phenopowerlaw
xi_0_sl: [31e6] xi_0_sl: [31e6]
xi_inf_sl: [63e6] xi_inf_sl: [63e6]
type: phenopowerlaw

23
processing/pre/geom_fromMinimalSurface.py
View File

@ -4,8 +4,6 @@ import os
import sys import sys
from optparse import OptionParser from optparse import OptionParser
import numpy as np
import damask import damask
@ -13,14 +11,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
minimal_surfaces = ['primitive','gyroid','diamond'] minimal_surfaces = list(damask.Geom._minimal_surface.keys())
surface = {
'primitive': lambda x,y,z: np.cos(x)+np.cos(y)+np.cos(z),
'gyroid': lambda x,y,z: np.sin(x)*np.cos(y)+np.sin(y)*np.cos(z)+np.cos(x)*np.sin(z),
'diamond': lambda x,y,z: np.cos(x-y)*np.cos(z)+np.sin(x+y)*np.sin(z),
}
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# MAIN # MAIN
@ -71,16 +62,8 @@ parser.set_defaults(type = minimal_surfaces[0],
name = None if filename == [] else filename[0] name = None if filename == [] else filename[0]
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
x,y,z = np.meshgrid(options.periods*2.0*np.pi*(np.arange(options.grid[0])+0.5)/options.grid[0], geom=damask.Geom.from_minimal_surface(options.grid,options.size,options.type,options.threshold,
options.periods*2.0*np.pi*(np.arange(options.grid[1])+0.5)/options.grid[1], options.periods,options.microstructure)
options.periods*2.0*np.pi*(np.arange(options.grid[2])+0.5)/options.grid[2],
indexing='xy',sparse=True)
microstructure = np.where(options.threshold < surface[options.type](x,y,z),
options.microstructure[1],options.microstructure[0])
geom=damask.Geom(microstructure,options.size,
comments=[scriptID + ' ' + ' '.join(sys.argv[1:])])
damask.util.croak(geom) damask.util.croak(geom)
geom.save_ASCII(sys.stdout if name is None else name,compress=False) geom.save_ASCII(sys.stdout if name is None else name,compress=False)

191
python/damask/_geom.py
View File

@ -120,6 +120,29 @@ class Geom:
return np.unique(self.material).size return np.unique(self.material).size
@staticmethod
def load(fname):
"""
Read a VTK rectilinear grid.
Parameters
----------
fname : str or or pathlib.Path
Geometry file to read.
Valid extension is .vtr, it will be appended if not given.
"""
v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
comments = v.get_comments()
grid = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
return Geom(material = v.get('material').reshape(grid,order='F'),
size = bbox[1] - bbox[0],
origin = bbox[0],
comments=comments)
@staticmethod @staticmethod
def load_ASCII(fname): def load_ASCII(fname):
""" """
@ -184,29 +207,6 @@ class Geom:
return Geom(material.reshape(grid,order='F'),size,origin,comments) return Geom(material.reshape(grid,order='F'),size,origin,comments)
@staticmethod
def load(fname):
"""
Read a VTK rectilinear grid.
Parameters
----------
fname : str or or pathlib.Path
Geometry file to read.
Valid extension is .vtr, it will be appended if not given.
"""
v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
comments = v.get_comments()
grid = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
return Geom(material = v.get('material').reshape(grid,order='F'),
size = bbox[1] - bbox[0],
origin = bbox[0],
comments=comments)
@staticmethod @staticmethod
def _find_closest_seed(seeds, weights, point): def _find_closest_seed(seeds, weights, point):
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights) return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
@ -292,6 +292,131 @@ class Geom:
) )
_minimal_surface = \
{'Schwarz P': lambda x,y,z: np.cos(x) + np.cos(y) + np.cos(z),
'Double Primitive': lambda x,y,z: ( 0.5 * (np.cos(x)*np.cos(y) + np.cos(y)*np.cos(z) + np.cos(z)*np.cos(x))
+ 0.2 * (np.cos(2*x) + np.cos(2*y) + np.cos(2*z)) ),
'Schwarz D': lambda x,y,z: ( np.sin(x)*np.sin(y)*np.sin(z)
+ np.sin(x)*np.cos(y)*np.cos(z)
+ np.cos(x)*np.cos(y)*np.sin(z)
+ np.cos(x)*np.sin(y)*np.cos(z) ),
'Complementary D': lambda x,y,z: ( np.cos(3*x+y)*np.cos(z) - np.sin(3*x-y)*np.sin(z) + np.cos(x+3*y)*np.cos(z)
+ np.sin(x-3*y)*np.sin(z) + np.cos(x-y)*np.cos(3*z) - np.sin(x+y)*np.sin(3*z) ),
'Double Diamond': lambda x,y,z: 0.5 * (np.sin(x)*np.sin(y)
+ np.sin(y)*np.sin(z)
+ np.sin(z)*np.sin(x)
+ np.cos(x) * np.cos(y) * np.cos(z) ),
'Dprime': lambda x,y,z: 0.5 * ( np.cos(x)*np.cos(y)*np.cos(z)
+ np.cos(x)*np.sin(y)*np.sin(z)
+ np.sin(x)*np.cos(y)*np.sin(z)
+ np.sin(x)*np.sin(y)*np.cos(z)
- np.sin(2*x)*np.sin(2*y)
- np.sin(2*y)*np.sin(2*z)
- np.sin(2*z)*np.sin(2*x) ) - 0.2,
'Gyroid': lambda x,y,z: np.cos(x)*np.sin(y) + np.cos(y)*np.sin(z) + np.cos(z)*np.sin(x),
'Gprime': lambda x,y,z : ( np.sin(2*x)*np.cos(y)*np.sin(z)
+ np.sin(2*y)*np.cos(z)*np.sin(x)
+ np.sin(2*z)*np.cos(x)*np.sin(y) ) + 0.32,
'Karcher K': lambda x,y,z: ( 0.3 * ( np.cos(x) + np.cos(y) + np.cos(z)
+ np.cos(x)*np.cos(y) + np.cos(y)*np.cos(z) + np.cos(z)*np.cos(x) )
- 0.4 * ( np.cos(2*x) + np.cos(2*y) + np.cos(2*z) ) ) + 0.2,
'Lidinoid': lambda x,y,z: 0.5 * ( np.sin(2*x)*np.cos(y)*np.sin(z)
+ np.sin(2*y)*np.cos(z)*np.sin(x)
+ np.sin(2*z)*np.cos(x)*np.sin(y)
- np.cos(2*x)*np.cos(2*y)
- np.cos(2*y)*np.cos(2*z)
- np.cos(2*z)*np.cos(2*x) ) + 0.15,
'Neovius': lambda x,y,z: ( 3 * (np.cos(x)+np.cos(y)+np.cos(z))
+ 4 * np.cos(x)*np.cos(y)*np.cos(z) ),
'Fisher-Koch S': lambda x,y,z: ( np.cos(2*x)*np.sin( y)*np.cos( z)
+ np.cos( x)*np.cos(2*y)*np.sin( z)
+ np.sin( x)*np.cos( y)*np.cos(2*z) ),
}
@staticmethod
def from_minimal_surface(grid,size,surface,threshold=0.0,periods=1,materials=(1,2)):
"""
Generate geometry from definition of triply periodic minimal surface.
Parameters
----------
grid : int numpy.ndarray of shape (3)
Number of grid points in x,y,z direction.
size : list or numpy.ndarray of shape (3)
Physical size of the geometry in meter.
surface : str
Type of the minimal surface. See notes for details.
threshold : float, optional.
Threshold of the minimal surface. Defaults to 0.0.
periods : integer, optional.
Number of periods per unit cell. Defaults to 1.
materials : (int, int), optional
Material IDs. Defaults to (1,2).
Notes
-----
The following triply-periodic minimal surfaces are implemented:
- Schwarz P
- Double Primitive
- Schwarz D
- Complementary D
- Double Diamond
- Dprime
- Gyroid
- Gprime
- Karcher K
- Lidinoid
- Neovius
- Fisher-Koch S
References
----------
Surface curvature in triply-periodic minimal surface architectures as
a distinct design parameter in preparing advanced tissue engineering scaffolds
Sébastien B G Blanquer, Maike Werner, Markus Hannula, Shahriar Sharifi,
Guillaume P R Lajoinie, David Eglin, Jari Hyttinen, André A Poot, and Dirk W Grijpma
10.1088/1758-5090/aa6553
Triply Periodic Bicontinuous Cubic Microdomain Morphologies by Symmetries
Meinhard Wohlgemuth, Nataliya Yufa, James Hoffman, and Edwin L. Thomas
10.1021/ma0019499
Minisurf – A minimal surface generator for finite element modeling and additive manufacturing
Meng-Ting Hsieh, Lorenzo Valdevit
10.1016/j.simpa.2020.100026
"""
x,y,z = np.meshgrid(periods*2.0*np.pi*(np.arange(grid[0])+0.5)/grid[0],
periods*2.0*np.pi*(np.arange(grid[1])+0.5)/grid[1],
periods*2.0*np.pi*(np.arange(grid[2])+0.5)/grid[2],
indexing='ij',sparse=True)
return Geom(material = np.where(threshold < Geom._minimal_surface[surface](x,y,z),materials[1],materials[0]),
size = size,
comments = util.execution_stamp('Geom','from_minimal_surface'),
)
def save(self,fname,compress=True):
"""
Generates vtk rectilinear grid.
Parameters
----------
fname : str, optional
Filename to write. If no file is given, a string is returned.
Valid extension is .vtr, it will be appended if not given.
compress : bool, optional
Compress with zlib algorithm. Defaults to True.
"""
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
v.add(self.material.flatten(order='F'),'material')
v.add_comments(self.comments)
v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress)
def save_ASCII(self,fname,compress=None): def save_ASCII(self,fname,compress=None):
""" """
Writes a geom file. Writes a geom file.
@ -364,26 +489,6 @@ class Geom:
f.write(f'{reps} of {former}\n') f.write(f'{reps} of {former}\n')
def save(self,fname,compress=True):
"""
Generates vtk rectilinear grid.
Parameters
----------
fname : str, optional
Filename to write. If no file is given, a string is returned.
Valid extension is .vtr, it will be appended if not given.
compress : bool, optional
Compress with zlib algorithm. Defaults to True.
"""
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
v.add(self.material.flatten(order='F'),'material')
v.add_comments(self.comments)
v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress)
def show(self): def show(self):
"""Show on screen.""" """Show on screen."""
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)

60
python/damask/_rotation.py
View File

@ -1,6 +1,8 @@
import numpy as np import numpy as np
from . import mechanics from . import mechanics
from . import util
from . import grid_filters
_P = -1 _P = -1
@ -647,13 +649,63 @@ class Rotation:
return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize() return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize()
# for compatibility (old names do not follow convention) # for compatibility
fromEulers = from_Eulers
fromQuaternion = from_quaternion
asAxisAngle = as_axis_angle
__mul__ = __matmul__ __mul__ = __matmul__
@staticmethod
def from_ODF(weights,Eulers,N=500,degrees=True,fractions=True,seed=None):
"""
Sample discrete values from a binned ODF.
Parameters
----------
weights : numpy.ndarray of shape (n)
Texture intensity values (probability density or volume fraction) at Euler grid points.
Eulers : numpy.ndarray of shape (n,3)
Grid coordinates in Euler space at which weights are defined.
N : integer, optional
Number of discrete orientations to be sampled from the given ODF.
Defaults to 500.
degrees : boolean, optional
Euler grid values are in degrees. Defaults to True.
fractions : boolean, optional
ODF values correspond to volume fractions, not probability density.
Defaults to True.
seed: {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
A seed to initialize the BitGenerator. Defaults to None, i.e. unpredictable entropy
will be pulled from the OS.
Returns
-------
samples : damask.Rotation of shape (N)
Array of sampled rotations closely representing the input ODF.
Notes
-----
Due to the distortion of Euler space in the vicinity of ϕ = 0, probability densities, p, defined on
grid points with ϕ = 0 will never result in reconstructed orientations as their dV/V = p dγ = p × 0.
Hence, it is recommended to transform any such dataset to cell centers that avoid grid points at ϕ = 0.
References
----------
P. Eisenlohr, F. Roters, Computational Materials Science 42(4), 670-678, 2008
https://doi.org/10.1016/j.commatsci.2007.09.015
"""
def _dg(eu,deg):
"""Return infinitesimal Euler space volume of bin(s)."""
Eulers_sorted = eu[np.lexsort((eu[:,0],eu[:,1],eu[:,2]))]
steps,size,_ = grid_filters.cell_coord0_gridSizeOrigin(Eulers_sorted)
delta = np.radians(size/steps) if deg else size/steps
return delta[0]*2.0*np.sin(delta[1]/2.0)*delta[2] / 8.0 / np.pi**2 * np.sin(np.radians(eu[:,1]) if deg else eu[:,1])
dg = 1.0 if fractions else _dg(Eulers,degrees)
dV_V = dg * np.maximum(0.0,weights.squeeze())
return Rotation.from_Eulers(Eulers[util.hybrid_IA(dV_V,N,seed)],degrees)
@staticmethod @staticmethod
def from_spherical_component(center,sigma,N=500,degrees=True,seed=None): def from_spherical_component(center,sigma,N=500,degrees=True,seed=None):
""" """

15
python/damask/util.py
View File

@ -20,6 +20,7 @@ __all__=[
'execute', 'execute',
'show_progress', 'show_progress',
'scale_to_coprime', 'scale_to_coprime',
'hybrid_IA',
'return_message', 'return_message',
'extendableOption', 'extendableOption',
'execution_stamp' 'execution_stamp'
@ -187,6 +188,20 @@ def execution_stamp(class_name,function_name=None):
return f'damask.{class_name}{_function_name} v{version} ({now})' return f'damask.{class_name}{_function_name} v{version} ({now})'
def hybrid_IA(dist,N,seed=None):
N_opt_samples,N_inv_samples = (max(np.count_nonzero(dist),N),0) # random subsampling if too little samples requested
scale_,scale,inc_factor = (0.0,float(N_opt_samples),1.0)
while (not np.isclose(scale, scale_)) and (N_inv_samples != N_opt_samples):
repeats = np.rint(scale*dist).astype(int)
N_inv_samples = np.sum(repeats)
scale_,scale,inc_factor = (scale,scale+inc_factor*0.5*(scale - scale_), inc_factor*2.0) \
if N_inv_samples < N_opt_samples else \
(scale_,0.5*(scale_ + scale), 1.0)
return np.repeat(np.arange(len(dist)),repeats)[np.random.default_rng(seed).permutation(N_inv_samples)[:N]]
#################################################################################################### ####################################################################################################
# Classes # Classes
#################################################################################################### ####################################################################################################

1524
python/tests/reference/Rotation/ODF_experimental.odf Normal file
View File

File diff suppressed because it is too large Load Diff

198510
python/tests/reference/Rotation/ODF_experimental.txt Normal file
View File

File diff suppressed because it is too large Load Diff

186629
python/tests/reference/Rotation/ODF_experimental_cell.txt Normal file
View File

File diff suppressed because it is too large Load Diff

18
python/tests/test_Colormap.py
View File

@ -77,17 +77,17 @@ class TestColormap:
@pytest.mark.parametrize('format',['ASCII','paraview','GOM','gmsh']) @pytest.mark.parametrize('format',['ASCII','paraview','GOM','gmsh'])
@pytest.mark.parametrize('model',['rgb','hsv','hsl','xyz','lab','msh']) @pytest.mark.parametrize('model',['rgb','hsv','hsl','xyz','lab','msh'])
def test_from_range(self,model,format,tmpdir): def test_from_range(self,model,format,tmp_path):
N = np.random.randint(2,256) N = np.random.randint(2,256)
c = Colormap.from_range(np.random.rand(3),np.random.rand(3),model=model,N=N) # noqa c = Colormap.from_range(np.random.rand(3),np.random.rand(3),model=model,N=N) # noqa
eval(f'c.save_{format}(tmpdir/"color_out")') eval(f'c.save_{format}(tmp_path/"color_out")')
@pytest.mark.parametrize('format',['ASCII','paraview','GOM','gmsh']) @pytest.mark.parametrize('format',['ASCII','paraview','GOM','gmsh'])
@pytest.mark.parametrize('name',['strain','gnuplot','Greys','PRGn','viridis']) @pytest.mark.parametrize('name',['strain','gnuplot','Greys','PRGn','viridis'])
def test_from_predefined(self,name,format,tmpdir): def test_from_predefined(self,name,format,tmp_path):
N = np.random.randint(2,256) N = np.random.randint(2,256)
c = Colormap.from_predefined(name,N) # noqa c = Colormap.from_predefined(name,N) # noqa
os.chdir(tmpdir) os.chdir(tmp_path)
eval(f'c.save_{format}()') eval(f'c.save_{format}()')
@pytest.mark.parametrize('format,name',[('ASCII','test.txt'), @pytest.mark.parametrize('format,name',[('ASCII','test.txt'),
@ -95,9 +95,9 @@ class TestColormap:
('GOM','test.legend'), ('GOM','test.legend'),
('gmsh','test.msh') ('gmsh','test.msh')
]) ])
def test_write_filehandle(self,format,name,tmpdir): def test_write_filehandle(self,format,name,tmp_path):
c = Colormap.from_predefined('Dark2') # noqa c = Colormap.from_predefined('Dark2') # noqa
fname = tmpdir/name fname = tmp_path/name
with open(fname,'w') as f: # noqa with open(fname,'w') as f: # noqa
eval(f'c.save_{format}(f)') eval(f'c.save_{format}(f)')
for i in range(10): for i in range(10):
@ -146,14 +146,14 @@ class TestColormap:
('GOM','.legend'), ('GOM','.legend'),
('gmsh','.msh') ('gmsh','.msh')
]) ])
def test_compare_reference(self,format,ext,tmpdir,reference_dir,update): def test_compare_reference(self,format,ext,tmp_path,reference_dir,update):
name = 'binary' name = 'binary'
c = Colormap.from_predefined(name) # noqa c = Colormap.from_predefined(name) # noqa
if update: if update:
os.chdir(reference_dir) os.chdir(reference_dir)
eval(f'c.save_{format}()') eval(f'c.save_{format}()')
else: else:
os.chdir(tmpdir) os.chdir(tmp_path)
eval(f'c.save_{format}()') eval(f'c.save_{format}()')
time.sleep(.5) time.sleep(.5)
assert filecmp.cmp(tmpdir/(name+ext),reference_dir/(name+ext)) assert filecmp.cmp(tmp_path/(name+ext),reference_dir/(name+ext))

82
python/tests/test_Geom.py
View File

@ -42,46 +42,46 @@ class TestGeom:
assert str(default.diff(new)) != '' assert str(default.diff(new)) != ''
def test_write_read_str(self,default,tmpdir): def test_write_read_str(self,default,tmp_path):
default.save_ASCII(str(tmpdir/'default.geom')) default.save_ASCII(str(tmp_path/'default.geom'))
new = Geom.load_ASCII(str(tmpdir/'default.geom')) new = Geom.load_ASCII(str(tmp_path/'default.geom'))
assert geom_equal(default,new) assert geom_equal(default,new)
def test_write_read_file(self,default,tmpdir): def test_write_read_file(self,default,tmp_path):
with open(tmpdir/'default.geom','w') as f: with open(tmp_path/'default.geom','w') as f:
default.save_ASCII(f,compress=True) default.save_ASCII(f,compress=True)
with open(tmpdir/'default.geom') as f: with open(tmp_path/'default.geom') as f:
new = Geom.load_ASCII(f) new = Geom.load_ASCII(f)
assert geom_equal(default,new) assert geom_equal(default,new)
def test_read_write_vtr(self,default,tmpdir): def test_read_write_vtr(self,default,tmp_path):
default.save(tmpdir/'default') default.save(tmp_path/'default')
for _ in range(10): for _ in range(10):
time.sleep(.2) time.sleep(.2)
if os.path.exists(tmpdir/'default.vtr'): break if os.path.exists(tmp_path/'default.vtr'): break
new = Geom.load(tmpdir/'default.vtr') new = Geom.load(tmp_path/'default.vtr')
assert geom_equal(new,default) assert geom_equal(new,default)
def test_invalid_geom(self,tmpdir): def test_invalid_geom(self,tmp_path):
with open('invalid_file','w') as f: with open(tmp_path/'invalid_file','w') as f:
f.write('this is not a valid header') f.write('this is not a valid header')
with open('invalid_file','r') as f: with open(tmp_path/'invalid_file','r') as f:
with pytest.raises(TypeError): with pytest.raises(TypeError):
Geom.load_ASCII(f) Geom.load_ASCII(f)
def test_invalid_vtr(self,tmpdir): def test_invalid_vtr(self,tmp_path):
v = VTK.from_rectilinearGrid(np.random.randint(5,10,3)*2,np.random.random(3) + 1.0) v = VTK.from_rectilinearGrid(np.random.randint(5,10,3)*2,np.random.random(3) + 1.0)
v.save(tmpdir/'no_materialpoint.vtr') v.save(tmp_path/'no_materialpoint.vtr')
for _ in range(10): for _ in range(10):
time.sleep(.2) time.sleep(.2)
if os.path.exists(tmpdir/'no_materialpoint.vtr'): break if os.path.exists(tmp_path/'no_materialpoint.vtr'): break
with pytest.raises(ValueError): with pytest.raises(ValueError):
Geom.load(tmpdir/'no_materialpoint.vtr') Geom.load(tmp_path/'no_materialpoint.vtr')
def test_invalid_material(self): def test_invalid_material(self):
with pytest.raises(TypeError): with pytest.raises(TypeError):
@ -92,9 +92,9 @@ class TestGeom:
assert g.material.dtype in np.sctypes['int'] assert g.material.dtype in np.sctypes['int']
@pytest.mark.parametrize('compress',[True,False]) @pytest.mark.parametrize('compress',[True,False])
def test_compress(self,default,tmpdir,compress): def test_compress(self,default,tmp_path,compress):
default.save_ASCII(tmpdir/'default.geom',compress=compress) default.save_ASCII(tmp_path/'default.geom',compress=compress)
new = Geom.load_ASCII(tmpdir/'default.geom') new = Geom.load_ASCII(tmp_path/'default.geom')
assert geom_equal(new,default) assert geom_equal(new,default)
@ -360,3 +360,45 @@ class TestGeom:
elif approach == 'Voronoi': elif approach == 'Voronoi':
geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5) geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5)
assert np.all(geom.material == material) assert np.all(geom.material == material)
@pytest.mark.parametrize('surface',['Schwarz P',
'Double Primitive',
'Schwarz D',
'Complementary D',
'Double Diamond',
'Dprime',
'Gyroid',
'Gprime',
'Karcher K',
'Lidinoid',
'Neovius',
'Fisher-Koch S',
])
def test_minimal_surface_basic_properties(self,surface):
grid = np.random.randint(60,100,3)
size = np.ones(3)+np.random.rand(3)
threshold = 2*np.random.rand()-1.
periods = np.random.randint(2)+1
materials = np.random.randint(0,40,2)
geom = Geom.from_minimal_surface(grid,size,surface,threshold,periods,materials)
assert set(geom.material.flatten()) | set(materials) == set(materials) \
and (geom.size == size).all() and (geom.grid == grid).all()
@pytest.mark.parametrize('surface,threshold',[('Schwarz P',0),
('Double Primitive',-1./6.),
('Schwarz D',0),
('Complementary D',0),
('Double Diamond',-0.133),
('Dprime',-0.0395),
('Gyroid',0),
('Gprime',0.22913),
('Karcher K',0.17045),
('Lidinoid',0.14455),
('Neovius',0),
('Fisher-Koch S',0),
])
def test_minimal_surface_volume(self,surface,threshold):
grid = np.ones(3,dtype=int)*64
geom = Geom.from_minimal_surface(grid,np.ones(3),surface,threshold)
assert np.isclose(np.count_nonzero(geom.material==1)/np.prod(geom.grid),.5,rtol=1e-3)

42
python/tests/test_Rotation.py
View File

@ -1,11 +1,11 @@
import os
import pytest import pytest
import numpy as np import numpy as np
from scipy import stats from scipy import stats
from damask import Rotation from damask import Rotation
from damask import Table
from damask import _rotation from damask import _rotation
from damask import grid_filters
n = 1000 n = 1000
atol=1.e-4 atol=1.e-4
@ -13,7 +13,7 @@ atol=1.e-4
@pytest.fixture @pytest.fixture
def reference_dir(reference_dir_base): def reference_dir(reference_dir_base):
"""Directory containing reference results.""" """Directory containing reference results."""
return os.path.join(reference_dir_base,'Rotation') return reference_dir_base/'Rotation'
@pytest.fixture @pytest.fixture
def set_of_rotations(set_of_quaternions): def set_of_rotations(set_of_quaternions):
@ -943,3 +943,39 @@ class TestRotation:
sigma_out = np.degrees(np.std(dist)) sigma_out = np.degrees(np.std(dist))
p = np.average(p) p = np.average(p)
assert (.9 < sigma/sigma_out < 1.1) and p > 1e-2, f'{sigma/sigma_out},{p}' assert (.9 < sigma/sigma_out < 1.1) and p > 1e-2, f'{sigma/sigma_out},{p}'
@pytest.mark.parametrize('fractions',[True,False])
@pytest.mark.parametrize('degrees',[True,False])
@pytest.mark.parametrize('N',[2**13,2**14,2**15])
def test_ODF_cell(self,reference_dir,fractions,degrees,N):
steps = np.array([144,36,36])
limits = np.array([360.,90.,90.])
rng = tuple(zip(np.zeros(3),limits))
weights = Table.load(reference_dir/'ODF_experimental_cell.txt').get('intensity').flatten()
Eulers = grid_filters.cell_coord0(steps,limits)
Eulers = np.radians(Eulers) if not degrees else Eulers
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees,fractions).as_Eulers(True)
weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
if fractions: assert np.sqrt(((weights_r - weights) ** 2).mean()) < 4
@pytest.mark.parametrize('degrees',[True,False])
@pytest.mark.parametrize('N',[2**13,2**14,2**15])
def test_ODF_node(self,reference_dir,degrees,N):
steps = np.array([144,36,36])
limits = np.array([360.,90.,90.])
rng = tuple(zip(-limits/steps*.5,limits-limits/steps*.5))
weights = Table.load(reference_dir/'ODF_experimental.txt').get('intensity')
weights = weights.reshape(steps+1,order='F')[:-1,:-1,:-1].reshape(-1,order='F')
Eulers = grid_filters.node_coord0(steps,limits)[:-1,:-1,:-1]
Eulers = np.radians(Eulers) if not degrees else Eulers
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees).as_Eulers(True)
weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5

24
python/tests/test_Table.py
View File

@ -34,31 +34,31 @@ class TestTable:
assert np.allclose(d,1.0) and d.shape[1:] == (1,) assert np.allclose(d,1.0) and d.shape[1:] == (1,)
@pytest.mark.parametrize('mode',['str','path']) @pytest.mark.parametrize('mode',['str','path'])
def test_write_read(self,default,tmpdir,mode): def test_write_read(self,default,tmp_path,mode):
default.save(tmpdir/'default.txt') default.save(tmp_path/'default.txt')
if mode == 'path': if mode == 'path':
new = Table.load(tmpdir/'default.txt') new = Table.load(tmp_path/'default.txt')
elif mode == 'str': elif mode == 'str':
new = Table.load(str(tmpdir/'default.txt')) new = Table.load(str(tmp_path/'default.txt'))
assert all(default.data==new.data) and default.shapes == new.shapes assert all(default.data==new.data) and default.shapes == new.shapes
def test_write_read_file(self,default,tmpdir): def test_write_read_file(self,default,tmp_path):
with open(tmpdir/'default.txt','w') as f: with open(tmp_path/'default.txt','w') as f:
default.save(f) default.save(f)
with open(tmpdir/'default.txt') as f: with open(tmp_path/'default.txt') as f:
new = Table.load(f) new = Table.load(f)
assert all(default.data==new.data) and default.shapes == new.shapes assert all(default.data==new.data) and default.shapes == new.shapes
def test_write_read_legacy_style(self,default,tmpdir): def test_write_read_legacy_style(self,default,tmp_path):
with open(tmpdir/'legacy.txt','w') as f: with open(tmp_path/'legacy.txt','w') as f:
default.save(f,legacy=True) default.save(f,legacy=True)
with open(tmpdir/'legacy.txt') as f: with open(tmp_path/'legacy.txt') as f:
new = Table.load(f) new = Table.load(f)
assert all(default.data==new.data) and default.shapes == new.shapes assert all(default.data==new.data) and default.shapes == new.shapes
def test_write_invalid_format(self,default,tmpdir): def test_write_invalid_format(self,default,tmp_path):
with pytest.raises(TypeError): with pytest.raises(TypeError):
default.save(tmpdir/'shouldnotbethere.txt',format='invalid') default.save(tmp_path/'shouldnotbethere.txt',format='invalid')
@pytest.mark.parametrize('mode',['str','path']) @pytest.mark.parametrize('mode',['str','path'])
def test_read_ang(self,reference_dir,mode): def test_read_ang(self,reference_dir,mode):

13
python/tests/test_util.py
View File

@ -1,5 +1,7 @@
import pytest import pytest
import numpy as np import numpy as np
from scipy import stats
from damask import util from damask import util
@ -31,3 +33,14 @@ class TestUtil:
def test_lackofprecision(self): def test_lackofprecision(self):
with pytest.raises(ValueError): with pytest.raises(ValueError):
util.scale_to_coprime(np.array([1/333.333,1,1])) util.scale_to_coprime(np.array([1/333.333,1,1]))
@pytest.mark.parametrize('rv',[stats.rayleigh(),stats.weibull_min(1.2),stats.halfnorm(),stats.pareto(2.62)])
def test_hybridIA(self,rv):
bins = np.linspace(0,10,100000)
centers = (bins[1:]+bins[:-1])/2
N_samples = bins.shape[0]-1000
dist = rv.pdf(centers)
selected = util.hybrid_IA(dist,N_samples)
dist_sampled = np.histogram(centers[selected],bins)[0]/N_samples*np.sum(dist)
assert np.sqrt(((dist - dist_sampled) ** 2).mean()) < .025 and selected.shape[0]==N_samples

4
src/CPFEM.f90
View File

@ -78,11 +78,11 @@ subroutine CPFEM_initAll
call DAMASK_interface_init call DAMASK_interface_init
call prec_init call prec_init
call IO_init call IO_init
call YAML_types_init
call YAML_parse_init
call config_init call config_init
call math_init call math_init
call rotations_init call rotations_init
call YAML_types_init
call YAML_parse_init
call HDF5_utilities_init call HDF5_utilities_init
call results_init(.false.) call results_init(.false.)
call discretization_marc_init call discretization_marc_init

4
src/CPFEM2.f90
View File

@ -48,11 +48,11 @@ subroutine CPFEM_initAll
#ifdef Mesh #ifdef Mesh
call FEM_quadrature_init call FEM_quadrature_init
#endif #endif
call YAML_types_init
call YAML_parse_init
call config_init call config_init
call math_init call math_init
call rotations_init call rotations_init
call YAML_types_init
call YAML_parse_init
call lattice_init call lattice_init
call HDF5_utilities_init call HDF5_utilities_init
call results_init(restart=interface_restartInc>0) call results_init(restart=interface_restartInc>0)

2
src/DAMASK_interface.f90
View File

@ -82,7 +82,7 @@ subroutine DAMASK_interface_init
print'(/,a)', ' <<<+- DAMASK_interface init -+>>>' print'(/,a)', ' <<<+- DAMASK_interface init -+>>>'
open(OUTPUT_unit, encoding='UTF-8') ! for special characters in output if(worldrank == 0) open(OUTPUT_UNIT, encoding='UTF-8') ! for special characters in output
! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK%203 ! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK%203
#ifdef DEBUG #ifdef DEBUG

17
src/IO.f90
View File

@ -427,20 +427,20 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
case (146) case (146)
msg = 'number of values does not match' msg = 'number of values does not match'
case (148) case (148)
msg = 'Nconstituents mismatch between homogenization and microstructure' msg = 'Nconstituents mismatch between homogenization and material'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! material error messages and related messages in mesh ! material error messages and related messages in mesh
case (150) case (150)
msg = 'index out of bounds' msg = 'index out of bounds'
case (151) case (151)
msg = 'microstructure has no constituents' msg = 'material has no constituents'
case (153) case (153)
msg = 'sum of phase fractions differs from 1' msg = 'sum of phase fractions differs from 1'
case (155) case (155)
msg = 'microstructure index out of bounds' msg = 'material index out of bounds'
case (180) case (180)
msg = 'missing/invalid microstructure definition via State Variable 2' msg = 'missing/invalid material definition via State Variable 2'
case (190) case (190)
msg = 'unknown element type:' msg = 'unknown element type:'
case (191) case (191)
@ -494,6 +494,10 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'Unsupported feature' msg = 'Unsupported feature'
case (706) case (706)
msg = 'Access by incorrect node type' msg = 'Access by incorrect node type'
case (707)
msg = 'Abrupt end of file'
case (708)
msg = '--- expected after YAML file header'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! errors related to the grid solver ! errors related to the grid solver
@ -522,7 +526,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
case (842) case (842)
msg = 'incomplete information in grid mesh header' msg = 'incomplete information in grid mesh header'
case (843) case (843)
msg = 'microstructure count mismatch' msg = 'material count mismatch'
case (844) case (844)
msg = 'invalid VTR file' msg = 'invalid VTR file'
case (846) case (846)
@ -621,6 +625,9 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
msg = 'polar decomposition failed' msg = 'polar decomposition failed'
case (700) case (700)
msg = 'unknown crystal symmetry' msg = 'unknown crystal symmetry'
case (709)
msg = 'read only the first document'
case (850) case (850)
msg = 'max number of cut back exceeded, terminating' msg = 'max number of cut back exceeded, terminating'
case default case default

120
src/YAML_parse.f90
View File

@ -227,6 +227,52 @@ logical function isKey(line)
end function isKey end function isKey
!--------------------------------------------------------------------------------------------------
! @brief skip empty lines
! @details update start position in the block by skipping empty lines if present.
!--------------------------------------------------------------------------------------------------
subroutine skip_empty_lines(blck,s_blck)
character(len=*), intent(in) :: blck
integer, intent(inout) :: s_blck
logical :: empty
empty = .true.
do while(empty .and. len_trim(blck(s_blck:)) /= 0)
empty = len_trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == 0
if(empty) s_blck = s_blck + index(blck(s_blck:),IO_EOL)
enddo
end subroutine skip_empty_lines
!--------------------------------------------------------------------------------------------------
! @brief skip file header
! @details update start position in the block by skipping file header if present.
!--------------------------------------------------------------------------------------------------
subroutine skip_file_header(blck,s_blck)
character(len=*), intent(in) :: blck
integer, intent(inout) :: s_blck
character(len=:), allocatable :: line
line = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))
if(index(adjustl(line),'%YAML') == 1) then
s_blck = s_blck + index(blck(s_blck:),IO_EOL)
call skip_empty_lines(blck,s_blck)
if(trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == '---') then
s_blck = s_blck + index(blck(s_blck:),IO_EOL)
else
call IO_error(708,ext_msg = line)
endif
endif
end subroutine skip_file_header
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! @brief reads a line of YAML block which is already in flow style ! @brief reads a line of YAML block which is already in flow style
! @details Dicts should be enlcosed within '{}' for it to be consistent with DAMASK YAML parser ! @details Dicts should be enlcosed within '{}' for it to be consistent with DAMASK YAML parser
@ -363,7 +409,9 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
do while (s_blck <= len_trim(blck)) do while (s_blck <= len_trim(blck))
e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2
line = IO_rmComment(blck(s_blck:e_blck)) line = IO_rmComment(blck(s_blck:e_blck))
if (len_trim(line) == 0) then if(trim(line) == '---' .or. trim(line) == '...') then
exit
elseif (len_trim(line) == 0) then
s_blck = e_blck + 2 ! forward to next line s_blck = e_blck + 2 ! forward to next line
cycle cycle
elseif(indentDepth(line,offset) > indent) then elseif(indentDepth(line,offset) > indent) then
@ -377,8 +425,10 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
else else
if(trim(adjustl(line)) == '-') then ! list item in next line if(trim(adjustl(line)) == '-') then ! list item in next line
s_blck = e_blck + 2 s_blck = e_blck + 2
e_blck = e_blck + index(blck(e_blck+2:),IO_EOL) call skip_empty_lines(blck,s_blck)
e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2
line = IO_rmComment(blck(s_blck:e_blck)) line = IO_rmComment(blck(s_blck:e_blck))
if(trim(line) == '---') call IO_error(707,ext_msg=line)
if(indentDepth(line) < indent .or. indentDepth(line) == indent) & if(indentDepth(line) < indent .or. indentDepth(line) == indent) &
call IO_error(701,ext_msg=line) call IO_error(701,ext_msg=line)
@ -447,7 +497,9 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset)
do while (s_blck <= len_trim(blck)) do while (s_blck <= len_trim(blck))
e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2
line = IO_rmComment(blck(s_blck:e_blck)) line = IO_rmComment(blck(s_blck:e_blck))
if (len_trim(line) == 0) then if(trim(line) == '---' .or. trim(line) == '...') then
exit
elseif (len_trim(line) == 0) then
s_blck = e_blck + 2 ! forward to next line s_blck = e_blck + 2 ! forward to next line
cycle cycle
elseif(indentDepth(line,offset) < indent) then elseif(indentDepth(line,offset) < indent) then
@ -510,10 +562,12 @@ recursive subroutine decide(blck,flow,s_blck,s_flow,offset)
character(len=:), allocatable :: line character(len=:), allocatable :: line
if(s_blck <= len(blck)) then if(s_blck <= len(blck)) then
call skip_empty_lines(blck,s_blck)
e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2
line = IO_rmComment(blck(s_blck:e_blck)) line = IO_rmComment(blck(s_blck:e_blck))
if(trim(line) == '---' .or. trim(line) == '...') then
if(len_trim(line) == 0) then continue ! end parsing at this point but not stop the simulation
elseif(len_trim(line) == 0) then
s_blck = e_blck +2 s_blck = e_blck +2
call decide(blck,flow,s_blck,s_flow,offset) call decide(blck,flow,s_blck,s_flow,offset)
elseif (isListItem(line)) then elseif (isListItem(line)) then
@ -548,23 +602,30 @@ function to_flow(blck)
character(len=:), allocatable :: to_flow character(len=:), allocatable :: to_flow
character(len=*), intent(in) :: blck !< YAML mixed style character(len=*), intent(in) :: blck !< YAML mixed style
character(len=:), allocatable :: line
integer :: s_blck, & !< start position in blck integer :: s_blck, & !< start position in blck
s_flow, & !< start position in flow s_flow, & !< start position in flow
offset, & !< counts leading '- ' in nested lists offset, & !< counts leading '- ' in nested lists
end_line end_line
if(isFlow(blck)) then
to_flow = trim(adjustl(blck)) allocate(character(len=len(blck)*2)::to_flow)
else s_flow = 1
allocate(character(len=len(blck)*2)::to_flow) s_blck = 1
! move forward here (skip empty lines) and remove '----' if found offset = 0
s_flow = 1
s_blck = 1 if(len_trim(blck) /= 0) then
offset = 0 call skip_empty_lines(blck,s_blck)
call skip_file_header(blck,s_blck)
line = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))
if(trim(line) == '---') s_blck = s_blck + index(blck(s_blck:),IO_EOL)
call decide(blck,to_flow,s_blck,s_flow,offset) call decide(blck,to_flow,s_blck,s_flow,offset)
to_flow = trim(to_flow(:s_flow-1))
endif endif
end_line = index(to_flow,IO_EOL) line = IO_rmComment(blck(s_blck:s_blck+index(blck(s_blck:),IO_EOL)-2))
if(end_line > 0) to_flow = to_flow(:end_line-1) if(trim(line)== '---') call IO_warning(709,ext_msg=line)
to_flow = trim(to_flow(:s_flow-1))
end_line = index(to_flow,IO_EOL)
if(end_line > 0) to_flow = to_flow(:end_line-1)
end function to_flow end function to_flow
@ -636,6 +697,20 @@ subroutine selfTest
if (.not. to_flow(block_dict_newline) == flow_dict) error stop 'to_flow' if (.not. to_flow(block_dict_newline) == flow_dict) error stop 'to_flow'
end block basic_dict end block basic_dict
only_flow: block
character(len=*), parameter :: flow_dict = &
" {a: [b,c: {d: e}, f: g, e]}"//IO_EOL
character(len=*), parameter :: flow_list = &
" [a,b: c, d,e: {f: g}]"//IO_EOL
character(len=*), parameter :: flow_1 = &
"{a: [b, {c: {d: e}}, {f: g}, e]}"
character(len=*), parameter :: flow_2 = &
"[a, {b: c}, d, {e: {f: g}}]"
if (.not. to_flow(flow_dict) == flow_1) error stop 'to_flow'
if (.not. to_flow(flow_list) == flow_2) error stop 'to_flow'
end block only_flow
basic_flow: block basic_flow: block
character(len=*), parameter :: flow_braces = & character(len=*), parameter :: flow_braces = &
" source: [{param: 1}, {param: 2}, {param: 3}, {param: 4}]"//IO_EOL " source: [{param: 1}, {param: 2}, {param: 3}, {param: 4}]"//IO_EOL
@ -650,12 +725,21 @@ subroutine selfTest
basic_mixed: block basic_mixed: block
character(len=*), parameter :: block_flow = & character(len=*), parameter :: block_flow = &
"%YAML 1.1"//IO_EOL//&
" "//IO_EOL//&
" "//IO_EOL//&
"---"//IO_EOL//&
" aa:"//IO_EOL//& " aa:"//IO_EOL//&
" - "//IO_EOL//& " - "//IO_EOL//&
" param_1: [a: b, c, {d: {e: [f: g, h]}}]"//IO_EOL//& " "//IO_EOL//&
" "//IO_EOL//&
" param_1: [a: b, c, {d: {e: [f: g, h]}}]"//IO_EOL//&
" - c: d"//IO_EOL//& " - c: d"//IO_EOL//&
" bb:"//IO_EOL//& " bb:"//IO_EOL//&
" - {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL " "//IO_EOL//&
" - "//IO_EOL//&
" {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL//&
"..."//IO_EOL
character(len=*), parameter :: mixed_flow = & character(len=*), parameter :: mixed_flow = &
"{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, {c: d}], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}" "{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, {c: d}], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}"

10
src/discretization.f90
View File

@ -15,7 +15,7 @@ module discretization
discretization_nElem discretization_nElem
integer, public, protected, dimension(:), allocatable :: & integer, public, protected, dimension(:), allocatable :: &
discretization_microstructureAt discretization_materialAt
real(pReal), public, protected, dimension(:,:), allocatable :: & real(pReal), public, protected, dimension(:,:), allocatable :: &
discretization_IPcoords0, & discretization_IPcoords0, &
@ -37,12 +37,12 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief stores the relevant information in globally accesible variables !> @brief stores the relevant information in globally accesible variables
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine discretization_init(microstructureAt,& subroutine discretization_init(materialAt,&
IPcoords0,NodeCoords0,& IPcoords0,NodeCoords0,&
sharedNodesBegin) sharedNodesBegin)
integer, dimension(:), intent(in) :: & integer, dimension(:), intent(in) :: &
microstructureAt materialAt
real(pReal), dimension(:,:), intent(in) :: & real(pReal), dimension(:,:), intent(in) :: &
IPcoords0, & IPcoords0, &
NodeCoords0 NodeCoords0
@ -51,10 +51,10 @@ subroutine discretization_init(microstructureAt,&
print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6) print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6)
discretization_nElem = size(microstructureAt,1) discretization_nElem = size(materialAt,1)
discretization_nIP = size(IPcoords0,2)/discretization_nElem discretization_nIP = size(IPcoords0,2)/discretization_nElem
discretization_microstructureAt = microstructureAt discretization_materialAt = materialAt
discretization_IPcoords0 = IPcoords0 discretization_IPcoords0 = IPcoords0
discretization_IPcoords = IPcoords0 discretization_IPcoords = IPcoords0

45
src/grid/DAMASK_grid.f90
View File

@ -181,7 +181,7 @@ program DAMASK_grid
if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check
call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase
newLoadCase%stress%myType='stress' newLoadCase%stress%myType='p'
field = 1 field = 1
newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default
thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
@ -210,18 +210,16 @@ program DAMASK_grid
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a * temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a *
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo enddo
newLoadCase%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) ! logical mask in 3x3 notation newLoadCase%deformation%mask = transpose(reshape(temp_maskVector,[ 3,3])) ! mask in 3x3 notation
newLoadCase%deformation%maskFloat = merge(ones,zeros,newLoadCase%deformation%maskLogical) ! float (1.0/0.0) mask in 3x3 notation newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation
newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation
case('p','stress', 's') case('p','stress', 's')
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
do j = 1, 9 do j = 1, 9
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo enddo
newLoadCase%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) newLoadCase%stress%mask = transpose(reshape(temp_maskVector,[ 3,3]))
newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical) newLoadCase%stress%values = math_9to33(temp_valueVector)
newLoadCase%stress%values = math_9to33(temp_valueVector)
case('t','time','delta') ! increment time case('t','time','delta') ! increment time
newLoadCase%time = IO_floatValue(line,chunkPos,i+1) newLoadCase%time = IO_floatValue(line,chunkPos,i+1)
case('n','incs','increments') ! number of increments case('n','incs','increments') ! number of increments
@ -268,8 +266,8 @@ program DAMASK_grid
print*, ' drop guessing along trajectory' print*, ' drop guessing along trajectory'
if (newLoadCase%deformation%myType == 'l') then if (newLoadCase%deformation%myType == 'l') then
do j = 1, 3 do j = 1, 3
if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & if (any(newLoadCase%deformation%mask(j,1:3) .eqv. .true.) .and. &
any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined any(newLoadCase%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined
enddo enddo
print*, ' velocity gradient:' print*, ' velocity gradient:'
else if (newLoadCase%deformation%myType == 'f') then else if (newLoadCase%deformation%myType == 'f') then
@ -278,20 +276,19 @@ program DAMASK_grid
print*, ' deformation gradient rate:' print*, ' deformation gradient rate:'
endif endif
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if(newLoadCase%deformation%maskLogical(i,j)) then if(newLoadCase%deformation%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j) write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j)
else else
write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' write(IO_STDOUT,'(2x,12a)',advance='no') ' * '
endif endif
enddo; write(IO_STDOUT,'(/)',advance='no') enddo; write(IO_STDOUT,'(/)',advance='no')
enddo enddo
if (any(newLoadCase%stress%maskLogical .eqv. & if (any(newLoadCase%stress%mask .eqv. newLoadCase%deformation%mask)) errorID = 831 ! exclusive or masking only
newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only if (any(newLoadCase%stress%mask .and. transpose(newLoadCase%stress%mask) .and. (math_I3<1))) &
if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) & errorID = 838 ! no rotation is allowed by stress BC
.and. (math_I3<1))) errorID = 838 ! no rotation is allowed by stress BC
print*, ' stress / GPa:' print*, ' stress / GPa:'
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if(newLoadCase%stress%maskLogical(i,j)) then if(newLoadCase%stress%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal
else else
write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' write(IO_STDOUT,'(2x,12a)',advance='no') ' * '
@ -368,11 +365,7 @@ program DAMASK_grid
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else else
if (currentLoadCase == 1) then ! 1st load case of logarithmic scale if (currentLoadCase == 1) then ! 1st load case of logarithmic scale
if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale timeinc = loadCases(1)%time*(2.0_pReal**real(max(inc-1,1)-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal))
endif
else ! not-1st load case of logarithmic scale else ! not-1st load case of logarithmic scale
timeinc = time0 * & timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/& ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/&
@ -432,17 +425,11 @@ program DAMASK_grid
do field = 1, nActiveFields do field = 1, nActiveFields
select case(loadCases(currentLoadCase)%ID(field)) select case(loadCases(currentLoadCase)%ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
solres(field) = mech_solution (& solres(field) = mech_solution(incInfo)
incInfo,timeinc,timeIncOld, &
stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rot)
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld) solres(field) = grid_thermal_spectral_solution(timeinc)
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld) solres(field) = grid_damage_spectral_solution(timeinc)
end select end select
if (.not. solres(field)%converged) exit ! no solution found if (.not. solres(field)%converged) exit ! no solution found

6
src/grid/grid_damage_spectral.f90
View File

@ -156,11 +156,10 @@ end subroutine grid_damage_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral damage scheme with internal iterations !> @brief solution for the spectral damage scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) function grid_damage_spectral_solution(timeinc) result(solution)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution timeinc !< increment in time for current solution
timeinc_old !< increment in time of last increment
integer :: i, j, k, cell integer :: i, j, k, cell
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
@ -174,7 +173,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)

126
src/grid/grid_mech_FEM.f90
View File

@ -25,63 +25,56 @@ module grid_mech_FEM
implicit none implicit none
private private
!-------------------------------------------------------------------------------------------------- type(tSolutionParams) :: params
! derived types
type(tSolutionParams), private :: params
type, private :: tNumerics type :: tNumerics
integer :: & integer :: &
itmin, & !< minimum number of iterations itmin, & !< minimum number of iterations
itmax !< maximum number of iterations itmax !< maximum number of iterations
real(pReal) :: & real(pReal) :: &
err_div, &
divTol, &
BCTol, &
eps_div_atol, & !< absolute tolerance for equilibrium eps_div_atol, & !< absolute tolerance for equilibrium
eps_div_rtol, & !< relative tolerance for equilibrium eps_div_rtol, & !< relative tolerance for equilibrium
eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC
eps_stress_rtol !< relative tolerance for fullfillment of stress BC eps_stress_rtol !< relative tolerance for fullfillment of stress BC
end type tNumerics end type tNumerics
type(tNumerics), private :: num type(tNumerics) :: num ! numerics parameters. Better name?
logical, private:: &
debugRotation logical :: debugRotation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM, private :: mech_grid DM :: mech_grid
SNES, private :: mech_snes SNES :: mech_snes
Vec, private :: solution_current, solution_lastInc, solution_rate Vec :: solution_current, solution_lastInc, solution_rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! common pointwise data ! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc real(pReal), dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc
real(pReal), private :: detJ real(pReal) :: detJ
real(pReal), private, dimension(3) :: delta real(pReal), dimension(3) :: delta
real(pReal), private, dimension(3,8) :: BMat real(pReal), dimension(3,8) :: BMat
real(pReal), private, dimension(8,8) :: HGMat real(pReal), dimension(8,8) :: HGMat
PetscInt, private :: xstart,ystart,zstart,xend,yend,zend PetscInt :: xstart,ystart,zstart,xend,yend,zend
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: & real(pReal), dimension(3,3) :: &
F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastIter = math_I3, &
F_aim_lastInc = math_I3, & !< previous average deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
P_aim = 0.0_pReal
character(len=:), allocatable, private :: incInfo !< time and increment information character(len=:), allocatable :: incInfo !< time and increment information
real(pReal), dimension(3,3,3,3) :: &
real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
S = 0.0_pReal !< current compliance (filled up with zeros) S = 0.0_pReal !< current compliance (filled up with zeros)
real(pReal), private :: & real(pReal) :: &
err_BC !< deviation from stress BC err_BC !< deviation from stress BC
integer, private :: & integer :: &
totalIter = 0 !< total iteration in current increment totalIter = 0 !< total iteration in current increment
public :: & public :: &
@ -99,7 +92,6 @@ contains
subroutine grid_mech_FEM_init subroutine grid_mech_FEM_init
real(pReal) :: HGCoeff = 0.0e-2_pReal real(pReal) :: HGCoeff = 0.0e-2_pReal
PetscInt, dimension(0:worldsize-1) :: localK
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal temp33_Real = 0.0_pReal
real(pReal), dimension(4,8) :: & real(pReal), dimension(4,8) :: &
@ -111,20 +103,21 @@ subroutine grid_mech_FEM_init
-1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, & -1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, &
1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, &
1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8])
real(pReal), dimension(3,3,3,3) :: devNull
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
PetscInt, dimension(0:worldsize-1) :: localK
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: & character(len=pStringLen) :: &
fileName fileName
class(tNode), pointer :: & class(tNode), pointer :: &
num_grid, & num_grid, &
debug_grid debug_grid
real(pReal), dimension(3,3,3,3) :: devNull
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT)
!----------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! debugging options ! debugging options
debug_grid => config_debug%get('grid', defaultVal=emptyList) debug_grid => config_debug%get('grid', defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation') debugRotation = debug_grid%contains('rotation')
@ -132,12 +125,12 @@ subroutine grid_mech_FEM_init
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict) num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal) num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) num%itmin = num_grid%get_asInt ('itmin', defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
@ -242,6 +235,7 @@ subroutine grid_mech_FEM_init
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity 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) F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restartRead endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(F) call utilities_updateCoords(F)
call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
@ -268,19 +262,12 @@ end subroutine grid_mech_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM scheme with internal iterations !> @brief solution for the FEM scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) function grid_mech_FEM_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
real(pReal), intent(in) :: &
timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
type(rotation), intent(in) :: &
rotation_BC
type(tSolutionState) :: & type(tSolutionState) :: &
solution solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -292,14 +279,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
@ -341,6 +321,14 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
PetscScalar, pointer, dimension(:,:,:,:) :: & PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc u_current,u_lastInc
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
@ -349,20 +337,20 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
else else
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = & F_aimDot = F_aimDot &
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = & F_aimDot = F_aimDot &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values + merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = & F_aimDot = F_aimDot &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
endif endif
if (guess) then if (guess) then
@ -382,6 +370,12 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc F_aim = F_aim_lastInc + F_aimDot * timeinc
if (stress_BC%myType=='p') then
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
elseif (stress_BC%myType=='pdot') then !UNTESTED
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
endif
call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr)
@ -497,8 +491,6 @@ subroutine formResidual(da_local,x_local, &
PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal
PetscScalar, dimension(8,3) :: x_elem, f_elem PetscScalar, dimension(8,3) :: x_elem, f_elem
PetscInt :: i, ii, j, jj, k, kk, ctr, ele PetscInt :: i, ii, j, jj, k, kk, ctr, ele
real(pReal), dimension(3,3) :: &
deltaF_aim
PetscInt :: & PetscInt :: &
PETScIter, & PETScIter, &
nfuncs nfuncs
@ -547,10 +539,8 @@ subroutine formResidual(da_local,x_local, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim_lastIter = F_aim F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask)))
F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual

100
src/grid/grid_mech_spectral_basic.f90
View File

@ -24,11 +24,9 @@ module grid_mech_spectral_basic
implicit none implicit none
private private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams) :: params type(tSolutionParams) :: params
type, private :: tNumerics type :: tNumerics
logical :: update_gamma !< update gamma operator with current stiffness logical :: update_gamma !< update gamma operator with current stiffness
integer :: & integer :: &
itmin, & !< minimum number of iterations itmin, & !< minimum number of iterations
@ -42,7 +40,7 @@ module grid_mech_spectral_basic
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
logical, private :: debugRotation logical :: debugRotation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
@ -62,10 +60,10 @@ module grid_mech_spectral_basic
F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
P_aim = 0.0_pReal
character(len=:), allocatable :: incInfo !< time and increment information character(len=:), allocatable :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
@ -96,18 +94,17 @@ subroutine grid_mech_spectral_basic_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal temp33_Real = 0.0_pReal
class (tNode), pointer :: &
num_grid, &
debug_grid
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: & PetscScalar, pointer, dimension(:,:,:,:) :: &
F ! pointer to solution data F ! pointer to solution data
PetscInt, dimension(worldsize) :: localK PetscInt, dimension(0:worldsize-1) :: localK
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
integer :: fileUnit integer :: fileUnit
character(len=pStringLen) :: & character(len=pStringLen) :: &
fileName fileName
class (tNode), pointer :: &
num_grid, &
debug_grid
print'(/,a)', ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(IO_STDOUT)
@ -126,13 +123,13 @@ subroutine grid_mech_spectral_basic_init
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict) num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.01_pReal) num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=1.0e-3_pReal)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) num%itmin = num_grid%get_asInt ('itmin',defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
@ -158,7 +155,7 @@ subroutine grid_mech_spectral_basic_init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
localK = 0 localK = 0
localK(worldrank+1) = grid3 localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
call DMDACreate3d(PETSC_COMM_WORLD, & call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
@ -202,8 +199,8 @@ subroutine grid_mech_spectral_basic_init
endif restartRead endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal) ! time increment 0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer
@ -231,19 +228,12 @@ end subroutine grid_mech_spectral_basic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the basic scheme with internal iterations !> @brief solution for the basic scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) function grid_mech_spectral_basic_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
real(pReal), intent(in) :: &
timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
type(rotation), intent(in) :: &
rotation_BC
type(tSolutionState) :: & type(tSolutionState) :: &
solution solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -255,17 +245,9 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
@ -303,7 +285,14 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
type(rotation), intent(in) :: & type(rotation), intent(in) :: &
rotation_BC rotation_BC
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F PetscScalar, pointer, dimension(:,:,:,:) :: F
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
@ -314,20 +303,20 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = & F_aimDot = F_aimDot &
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = & F_aimDot = F_aimDot &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values + merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = & F_aimDot = F_aimDot &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
endif endif
Fdot = utilities_calculateRate(guess, & Fdot = utilities_calculateRate(guess, &
@ -341,7 +330,13 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc F_aim = F_aim_lastInc + F_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average if (stress_BC%myType=='p') then
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
elseif (stress_BC%myType=='pdot') then !UNTESTED
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
endif
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
@ -375,7 +370,7 @@ subroutine grid_mech_spectral_basic_restartWrite
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
print'(a)', ' writing solver data required for restart to file'; flush(IO_STDOUT) print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'w') fileHandle = HDF5_openFile(fileName,'w')
@ -469,6 +464,7 @@ subroutine formResidual(in, F, &
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
@ -491,16 +487,16 @@ subroutine formResidual(in, F, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
F_aim = F_aim - deltaF_aim F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask)))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme ! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier
call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier

107
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -15,7 +15,6 @@ module grid_mech_spectral_polarisation
use DAMASK_interface use DAMASK_interface
use HDF5_utilities use HDF5_utilities
use math use math
use rotations
use spectral_utilities use spectral_utilities
use FEsolving use FEsolving
use config use config
@ -25,8 +24,6 @@ module grid_mech_spectral_polarisation
implicit none implicit none
private private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams) :: params type(tSolutionParams) :: params
type :: tNumerics type :: tNumerics
@ -48,7 +45,7 @@ module grid_mech_spectral_polarisation
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
logical, private :: debugRotation logical :: debugRotation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
@ -71,8 +68,8 @@ module grid_mech_spectral_polarisation
F_aim = math_I3, & !< current prescribed deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient
F_av = 0.0_pReal, & !< average incompatible def grad field F_av = 0.0_pReal, & !< average incompatible def grad field
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
P_aim = 0.0_pReal
character(len=:), allocatable :: incInfo !< time and increment information character(len=:), allocatable :: incInfo !< time and increment information
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvg = 0.0_pReal, & !< current volume average stiffness
@ -108,10 +105,6 @@ subroutine grid_mech_spectral_polarisation_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal temp33_Real = 0.0_pReal
class (tNode), pointer :: &
num_grid, &
debug_grid
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: & PetscScalar, pointer, dimension(:,:,:,:) :: &
FandF_tau, & ! overall pointer to solution data FandF_tau, & ! overall pointer to solution data
@ -122,13 +115,16 @@ subroutine grid_mech_spectral_polarisation_init
integer :: fileUnit integer :: fileUnit
character(len=pStringLen) :: & character(len=pStringLen) :: &
fileName fileName
class (tNode), pointer :: &
num_grid, &
debug_grid
print'(/,a)', ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(IO_STDOUT)
print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
!------------------------------------------------------------------------------------------------ !-------------------------------------------------------------------------------------------------
! debugging options ! debugging options
debug_grid => config_debug%get('grid',defaultVal=emptyList) debug_grid => config_debug%get('grid',defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation') debugRotation = debug_grid%contains('rotation')
@ -136,17 +132,18 @@ subroutine grid_mech_spectral_polarisation_init
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict) num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_curl_atol = num_grid%get_asFloat ('eps_curl_atol', defaultVal=1.0e-10_pReal) num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_curl_rtol = num_grid%get_asFloat ('eps_curl_rtol', defaultVal=5.0e-4_pReal) num%eps_curl_atol = num_grid%get_asFloat('eps_curl_atol', defaultVal=1.0e-10_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal) num%eps_curl_rtol = num_grid%get_asFloat('eps_curl_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) num%itmin = num_grid%get_asInt ('itmin', defaultVal=1)
num%alpha = num_grid%get_asFloat ('alpha', defaultVal=1.0_pReal) num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
num%beta = num_grid%get_asFloat ('beta', defaultVal=1.0_pReal) num%alpha = num_grid%get_asFloat('alpha', defaultVal=1.0_pReal)
num%beta = num_grid%get_asFloat('beta', defaultVal=1.0_pReal)
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
@ -228,8 +225,8 @@ subroutine grid_mech_spectral_polarisation_init
endif restartRead endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal) ! time increment 0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
@ -259,19 +256,12 @@ end subroutine grid_mech_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the Polarisation scheme with internal iterations !> @brief solution for the Polarisation scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
real(pReal), intent(in) :: &
timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
type(rotation), intent(in) :: &
rotation_BC
type(tSolutionState) :: & type(tSolutionState) :: &
solution solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -283,21 +273,13 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
if (num%update_gamma) then if(num%update_gamma) then
call utilities_updateGamma(C_minMaxAvg) call utilities_updateGamma(C_minMaxAvg)
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg) S_scale = math_invSym3333(C_minMaxAvg)
endif endif
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
@ -335,10 +317,17 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
type(rotation), intent(in) :: & type(rotation), intent(in) :: &
rotation_BC rotation_BC
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau
integer :: i, j, k integer :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33 real(pReal), dimension(3,3) :: F_lambda33
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
F => FandF_tau(0: 8,:,:,:) F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:) F_tau => FandF_tau(9:17,:,:,:)
@ -350,20 +339,20 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = & F_aimDot = F_aimDot &
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = & F_aimDot = F_aimDot &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values + merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = & F_aimDot = F_aimDot &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
endif endif
Fdot = utilities_calculateRate(guess, & Fdot = utilities_calculateRate(guess, &
@ -381,6 +370,12 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc F_aim = F_aim_lastInc + F_aimDot * timeinc
if (stress_BC%myType=='p') then
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
elseif (stress_BC%myType=='pdot') then !UNTESTED
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
endif
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),& rotation_BC%rotate(F_aim,active=.true.)),&
[9,grid(1),grid(2),grid3]) [9,grid(1),grid(2),grid3])
@ -595,15 +590,15 @@ subroutine formResidual(in, FandF_tau, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & err_BC = maxval(abs(merge(P_av-P_aim, &
-params%rotation_BC%rotate(F_av)) + & math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),&
params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc params%stress_mask)))
! calculate divergence ! calculate divergence
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise
call utilities_FFTtensorForward call utilities_FFTtensorForward
err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
@ -622,7 +617,7 @@ subroutine formResidual(in, FandF_tau, &
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward call utilities_FFTtensorForward
err_curl = Utilities_curlRMS() err_curl = utilities_curlRMS()
end subroutine formResidual end subroutine formResidual

6
src/grid/grid_thermal_spectral.f90
View File

@ -148,11 +148,10 @@ end subroutine grid_thermal_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral thermal scheme with internal iterations !> @brief solution for the spectral thermal scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) function grid_thermal_spectral_solution(timeinc) result(solution)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution timeinc !< increment in time for current solution
timeinc_old !< increment in time of last increment
integer :: i, j, k, cell integer :: i, j, k, cell
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
@ -166,7 +165,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old
call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr)

31
src/grid/spectral_utilities.f90
View File

@ -47,15 +47,15 @@ module spectral_utilities
complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field fourier representation for fftw complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field fourier representation for fftw
real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for fftw real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for fftw
complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field fourier representation for fftw complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field fourier representation for fftw
complex(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method complex(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method
complex(pReal), private, dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives complex(pReal), dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives
complex(pReal), private, dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives complex(pReal), dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives
real(pReal), private, dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness real(pReal), dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! plans for FFTW ! plans for FFTW
type(C_PTR), private :: & type(C_PTR) :: &
planTensorForth, & !< FFTW MPI plan P(x) to P(k) planTensorForth, & !< FFTW MPI plan P(x) to P(k)
planTensorBack, & !< FFTW MPI plan F(k) to F(x) planTensorBack, & !< FFTW MPI plan F(k) to F(x)
planVectorForth, & !< FFTW MPI plan v(x) to v(k) planVectorForth, & !< FFTW MPI plan v(x) to v(k)
@ -65,7 +65,7 @@ module spectral_utilities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables controlling debugging ! variables controlling debugging
logical, private :: & logical :: &
debugGeneral, & !< general debugging of spectral solver debugGeneral, & !< general debugging of spectral solver
debugRotation, & !< also printing out results in lab frame debugRotation, & !< also printing out results in lab frame
debugPETSc !< use some in debug defined options for more verbose PETSc solution debugPETSc !< use some in debug defined options for more verbose PETSc solution
@ -82,10 +82,9 @@ module spectral_utilities
end type tSolutionState end type tSolutionState
type, public :: tBoundaryCondition !< set of parameters defining a boundary condition type, public :: tBoundaryCondition !< set of parameters defining a boundary condition
real(pReal), dimension(3,3) :: values = 0.0_pReal, & real(pReal), dimension(3,3) :: values = 0.0_pReal
maskFloat = 0.0_pReal logical, dimension(3,3) :: mask = .false.
logical, dimension(3,3) :: maskLogical = .false. character(len=pStringLen) :: myType = 'None'
character(len=pStringLen) :: myType = 'None'
end type tBoundaryCondition end type tBoundaryCondition
type, public :: tLoadCase type, public :: tLoadCase
@ -101,21 +100,21 @@ module spectral_utilities
integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
end type tLoadCase end type tLoadCase
type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase type, public :: tSolutionParams
real(pReal), dimension(3,3) :: stress_mask, stress_BC real(pReal), dimension(3,3) :: stress_BC
logical, dimension(3,3) :: stress_mask
type(rotation) :: rotation_BC type(rotation) :: rotation_BC
real(pReal) :: timeinc real(pReal) :: timeinc, timeincOld
real(pReal) :: timeincOld
end type tSolutionParams end type tSolutionParams
type, private :: tNumerics type :: tNumerics
integer :: & integer :: &
divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints] divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints]
logical :: & logical :: &
memory_efficient !< calculate gamma operator on the fly memory_efficient !< calculate gamma operator on the fly
end type tNumerics end type tNumerics
type(tNumerics), private :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
DERIVATIVE_CONTINUOUS_ID, & DERIVATIVE_CONTINUOUS_ID, &

50
src/material.f90
View File

@ -2,7 +2,7 @@
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Parses material config file, either solverJobName.materialConfig or material.config !> @brief Defines phase and homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module material module material
use prec use prec
@ -97,7 +97,7 @@ module material
material_orientation0 !< initial orientation of each grain,IP,element material_orientation0 !< initial orientation of each grain,IP,element
integer, dimension(:), allocatable, private :: & integer, dimension(:), allocatable, private :: &
microstructure_Nconstituents !< number of constituents in each microstructure material_Nconstituents !< number of constituents in each material
@ -180,8 +180,8 @@ subroutine material_init(restart)
material_name_homogenization(myHomog) = trim(adjustl(sectionName))//material_homogenization%getKey(myHomog) material_name_homogenization(myHomog) = trim(adjustl(sectionName))//material_homogenization%getKey(myHomog)
enddo enddo
call material_parseMicrostructure call material_parseMaterial
print*, 'Microstructure parsed' print*, 'Material parsed'
call material_parseHomogenization call material_parseHomogenization
print*, 'Homogenization parsed' print*, 'Homogenization parsed'
@ -317,16 +317,16 @@ end subroutine material_parseHomogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief parses the microstructure part in the material configuration file !> @brief parses the material part in the material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine material_parseMicrostructure subroutine material_parseMaterial
class(tNode), pointer :: microstructures, & !> list of microstructures class(tNode), pointer :: materials, & !> list of materials
microstructure, & !> microstructure definition material, & !> material definition
constituents, & !> list of constituents constituents, & !> list of constituents
constituent, & !> constituent definition constituent, & !> constituent definition
phases, & phases, &
homogenization homogenizations
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
counterPhase, & counterPhase, &
@ -341,17 +341,17 @@ subroutine material_parseMicrostructure
c, & c, &
maxNconstituents maxNconstituents
microstructures => config_material%get('microstructure') materials => config_material%get('material')
if(any(discretization_microstructureAt > microstructures%length)) & if(any(discretization_materialAt > materials%length)) &
call IO_error(155,ext_msg='More microstructures requested than found in material.yaml') call IO_error(155,ext_msg='More materials requested than found in material.yaml')
allocate(microstructure_Nconstituents(microstructures%length),source=0) allocate(material_Nconstituents(materials%length),source=0)
do m = 1, microstructures%length do m = 1, materials%length
microstructure => microstructures%get(m) material => materials%get(m)
constituents => microstructure%get('constituents') constituents => material%get('constituents')
microstructure_Nconstituents(m) = constituents%length material_Nconstituents(m) = constituents%length
enddo enddo
maxNconstituents = maxval(microstructure_Nconstituents) maxNconstituents = maxval(material_Nconstituents)
allocate(material_homogenizationAt(discretization_nElem),source=0) allocate(material_homogenizationAt(discretization_nElem),source=0)
allocate(material_homogenizationMemberAt(discretization_nIP,discretization_nElem),source=0) allocate(material_homogenizationMemberAt(discretization_nIP,discretization_nElem),source=0)
@ -362,14 +362,14 @@ subroutine material_parseMicrostructure
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(counterPhase(phases%length),source=0) allocate(counterPhase(phases%length),source=0)
homogenization => config_material%get('homogenization') homogenizations => config_material%get('homogenization')
allocate(counterHomogenization(homogenization%length),source=0) allocate(counterHomogenization(homogenizations%length),source=0)
do e = 1, discretization_nElem do e = 1, discretization_nElem
microstructure => microstructures%get(discretization_microstructureAt(e)) material => materials%get(discretization_materialAt(e))
constituents => microstructure%get('constituents') constituents => material%get('constituents')
material_homogenizationAt(e) = homogenization%getIndex(microstructure%get_asString('homogenization')) material_homogenizationAt(e) = homogenizations%getIndex(material%get_asString('homogenization'))
do i = 1, discretization_nIP do i = 1, discretization_nIP
counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1 counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1
material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e)) material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e))
@ -385,7 +385,7 @@ subroutine material_parseMicrostructure
counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1 counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1
material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e)) material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e))
call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('orientation',requiredSize=4)) call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4))
enddo enddo
enddo enddo
@ -393,7 +393,7 @@ subroutine material_parseMicrostructure
enddo enddo
end subroutine material_parseMicrostructure end subroutine material_parseMaterial
end module material end module material

8
src/mesh/discretization_mesh.f90
View File

@ -78,7 +78,7 @@ subroutine discretization_mesh_init(restart)
IS :: faceSetIS IS :: faceSetIS
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
microstructureAt materialAt
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
integer :: integrationOrder !< order of quadrature rule required integer :: integrationOrder !< order of quadrature rule required
@ -164,9 +164,9 @@ subroutine discretization_mesh_init(restart)
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p)
call mesh_FEM_build_ipVolumes(dimPlex) call mesh_FEM_build_ipVolumes(dimPlex)
allocate(microstructureAt(mesh_NcpElems)) allocate(materialAt(mesh_NcpElems))
do j = 1, mesh_NcpElems do j = 1, mesh_NcpElems
call DMGetLabelValue(geomMesh,'material',j-1,microstructureAt(j),ierr) call DMGetLabelValue(geomMesh,'material',j-1,materialAt(j),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
end do end do
@ -178,7 +178,7 @@ subroutine discretization_mesh_init(restart)
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
call discretization_init(microstructureAt,& call discretization_init(materialAt,&
reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), &
mesh_node0) mesh_node0)

2
src/rotations.f90
View File

@ -1452,7 +1452,7 @@ subroutine selfTest
contains contains
function quaternion_equal(qu1,qu2) result(ok) pure recursive function quaternion_equal(qu1,qu2) result(ok)
real(pReal), intent(in), dimension(4) :: qu1,qu2 real(pReal), intent(in), dimension(4) :: qu1,qu2
logical :: ok logical :: ok