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 dislotwin-fix-TWIP-TRIP

This commit is contained in:
Martin Diehl 2022-02-11 06:05:09 +01:00
parent 8f0a16e654 5f7a18c32d
commit 90dabadf08
64 changed files with 1512 additions and 1220 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
.gitlab-ci.yml
View File

@ -36,21 +36,21 @@ variables:
# Names of module files to load
# ===============================================================================================
# ++++++++++++ Compiler +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
COMPILER_GNU: "Compiler/GNU/10"
COMPILER_GNU: "Compiler/GNU/10"
COMPILER_INTELLLVM: "Compiler/oneAPI/2022.0.1 Libraries/IMKL/2022.0.1"
COMPILER_INTEL: "Compiler/Intel/2022.0.1 Libraries/IMKL/2022.0.1"
# ++++++++++++ MPI ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.1"
MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0"
MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.2"
MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0"
MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
# ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
PETSC_GNU: "Libraries/PETSc/3.16.4/GNU-10-OpenMPI-4.1.2"
PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.3/Intel-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.4/Intel-2022.0.1-IntelMPI-2021.5.0"
# ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MSC: "FEM/MSC/2021.3.1"
IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
HDF5Marc: "HDF5/1.12.1/Intel-19.1.2"
MSC: "FEM/MSC/2021.3.1"
IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
HDF5Marc: "HDF5/1.12.1/Intel-19.1.2"
###################################################################################################

2
CMakeLists.txt
View File

@ -42,7 +42,7 @@ string(TOUPPER "${CMAKE_BUILD_TYPE}" CMAKE_BUILD_TYPE)
if(CMAKE_BUILD_TYPE STREQUAL "DEBUG" OR CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
set(DEBUG_FLAGS "${DEBUG_FLAGS} -DDEBUG")
set(PARALLEL "OFF")
set(OPTI "OFF")
set(OPTI "DEBUG")
elseif(CMAKE_BUILD_TYPE STREQUAL "RELEASE")
set(PARALLEL "ON")
set(OPTI "DEFENSIVE")

2
PRIVATE

@ -1 +1 @@
Subproject commit 5598ec4892b0921fccf63e8570f9fcd8e0365716
Subproject commit 5774122bf48d637704bb4afb10b87c34a4dbcaba

17
cmake/Compiler-GNU.cmake
View File

@ -9,26 +9,24 @@ if (OPENMP)
set (OPENMP_FLAGS "-fopenmp")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
if (OPTIMIZATION STREQUAL "DEBUG")
set (OPTIMIZATION_FLAGS "-Og")
elseif (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2 -mtune=generic")
set (OPTIMIZATION_FLAGS "-O2 -mtune=native")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-O3 -march=native -ffast-math -funroll-loops -ftree-vectorize")
set (OPTIMIZATION_FLAGS "-O3 -march=native -funroll-loops -ftree-vectorize -flto")
endif ()
set (STANDARD_CHECK "-std=f2018 -pedantic-errors" )
set (LINKER_FLAGS "${LINKER_FLAGS} -Wl")
# options parsed directly to the linker
set (LINKER_FLAGS "${LINKER_FLAGS},-undefined,dynamic_lookup" )
# ensure to link against dynamic libraries
#------------------------------------------------------------------------------------------------
# Fine tuning compilation options
set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp")
# preprocessor
set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIC -fPIE")
set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIE")
# position independent code
set (COMPILE_FLAGS "${COMPILE_FLAGS} -ffree-line-length-132")
@ -123,6 +121,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -ffpe-trap=invalid,zero,overflow")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -g")
# Generate symbolic debugging information in the object file
set (DEBUG_FLAGS "${DEBUG_FLAGS} -Og")
# Optimize debugging experience
set (DEBUG_FLAGS "${DEBUG_FLAGS} -fbacktrace")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -fdump-core")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -fcheck=all")

7
cmake/Compiler-Intel.cmake
View File

@ -9,12 +9,12 @@ if (OPENMP)
set (OPENMP_FLAGS "-qopenmp -parallel")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG")
set (OPTIMIZATION_FLAGS "-O0 -no-ip")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost")
set (OPTIMIZATION_FLAGS "-ipo -O3 -fp-model fast=2 -xHost")
# -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost"
endif ()
@ -110,6 +110,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0")
# generate debug information for parameters
# Disabled due to ICE when compiling phase_damage.f90 (not understandable, there is no parameter in there)
set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
# generate complete debugging information
# Additional options
# -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
# -check: Checks at runtime, where

5
cmake/Compiler-IntelLLVM.cmake
View File

@ -9,7 +9,7 @@ if (OPENMP)
set (OPENMP_FLAGS "-qopenmp")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG")
set (OPTIMIZATION_FLAGS "-O0")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
@ -109,6 +109,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all")
# generate debug information for parameters
set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
# generate complete debugging information
# Additional options
# -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
# -check: Checks at runtime, where

178
processing/legacy/geom_grainGrowth.py vendored
View File

@ -1,178 +0,0 @@
#!/usr/bin/env python3
import os
import sys
from io import StringIO
from optparse import OptionParser
import numpy as np
from scipy import ndimage
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
getInterfaceEnergy = lambda A,B: np.float32((A != B)*1.0) # 1.0 if A & B are distinct, 0.0 otherwise
struc = ndimage.generate_binary_structure(3,1) # 3D von Neumann neighborhood
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [geomfile(s)]', description = """
Smoothen interface roughness by simulated curvature flow.
This is achieved by the diffusion of each initially sharply bounded grain volume within the periodic domain
up to a given distance 'd' voxels.
The final geometry is assembled by selecting at each voxel that grain index for which the concentration remains largest.
References 10.1073/pnas.1111557108 (10.1006/jcph.1994.1105)
""", version = scriptID)
parser.add_option('-d', '--distance',
dest = 'd',
type = 'float', metavar = 'float',
help = 'diffusion distance in voxels [%default]')
parser.add_option('-N', '--iterations',
dest = 'N',
type = 'int', metavar = 'int',
help = 'curvature flow iterations [%default]')
parser.add_option('-i', '--immutable',
action = 'extend', dest = 'immutable', metavar = '<int LIST>',
help = 'list of immutable material indices')
parser.add_option('--ndimage',
dest = 'ndimage', action='store_true',
help = 'use ndimage.gaussian_filter in lieu of explicit FFT')
parser.set_defaults(d = 1,
N = 1,
immutable = [],
ndimage = False,
)
(options, filenames) = parser.parse_args()
options.immutable = list(map(int,options.immutable))
if filenames == []: filenames = [None]
for name in filenames:
damask.util.report(scriptName,name)
geom = damask.Grid.load(StringIO(''.join(sys.stdin.read())) if name is None else name)
grid_original = geom.cells
damask.util.croak(geom)
material = np.tile(geom.material,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1
grid = np.array(material.shape)
# --- initialize support data ---------------------------------------------------------------------
# store a copy of the initial material indices to find locations of immutable indices
material_original = np.copy(material)
if not options.ndimage:
X,Y,Z = np.mgrid[0:grid[0],0:grid[1],0:grid[2]]
# Calculates gaussian weights for simulating 3d diffusion
gauss = np.exp(-(X*X + Y*Y + Z*Z)/(2.0*options.d*options.d),dtype=np.float32) \
/np.power(2.0*np.pi*options.d*options.d,(3.0 - np.count_nonzero(grid_original == 1))/2.,dtype=np.float32)
gauss[:,:,:grid[2]//2:-1] = gauss[:,:,1:(grid[2]+1)//2] # trying to cope with uneven (odd) grid size
gauss[:,:grid[1]//2:-1,:] = gauss[:,1:(grid[1]+1)//2,:]
gauss[:grid[0]//2:-1,:,:] = gauss[1:(grid[0]+1)//2,:,:]
gauss = np.fft.rfftn(gauss).astype(np.complex64)
for smoothIter in range(options.N):
interfaceEnergy = np.zeros(material.shape,dtype=np.float32)
for i in (-1,0,1):
for j in (-1,0,1):
for k in (-1,0,1):
# assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood)
interfaceEnergy = np.maximum(interfaceEnergy,
getInterfaceEnergy(material,np.roll(np.roll(np.roll(
material,i,axis=0), j,axis=1), k,axis=2)))
# periodically extend interfacial energy array by half a grid size in positive and negative directions
periodic_interfaceEnergy = np.tile(interfaceEnergy,(3,3,3))[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2]
# transform bulk volume (i.e. where interfacial energy remained zero), store index of closest boundary voxel
index = ndimage.morphology.distance_transform_edt(periodic_interfaceEnergy == 0.,
return_distances = False,
return_indices = True)
# want array index of nearest voxel on periodically extended boundary
periodic_bulkEnergy = periodic_interfaceEnergy[index[0],
index[1],
index[2]].reshape(2*grid) # fill bulk with energy of nearest interface
if options.ndimage:
periodic_diffusedEnergy = ndimage.gaussian_filter(
np.where(ndimage.morphology.binary_dilation(periodic_interfaceEnergy > 0.,
structure = struc,
iterations = int(round(options.d*2.))-1, # fat boundary
),
periodic_bulkEnergy, # ...and zero everywhere else
0.),
sigma = options.d)
else:
diffusedEnergy = np.fft.irfftn(np.fft.rfftn(
np.where(
ndimage.morphology.binary_dilation(interfaceEnergy > 0.,
structure = struc,
iterations = int(round(options.d*2.))-1),# fat boundary
periodic_bulkEnergy[grid[0]//2:-grid[0]//2, # retain filled energy on fat boundary...
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2], # ...and zero everywhere else
0.)).astype(np.complex64) *
gauss).astype(np.float32)
periodic_diffusedEnergy = np.tile(diffusedEnergy,(3,3,3))[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # periodically extend the smoothed bulk energy
# transform voxels close to interface region
index = ndimage.morphology.distance_transform_edt(periodic_diffusedEnergy >= 0.95*np.amax(periodic_diffusedEnergy),
return_distances = False,
return_indices = True) # want index of closest bulk grain
periodic_material = np.tile(material,(3,3,3))[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # periodically extend the geometry
material = periodic_material[index[0],
index[1],
index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # extent grains into interface region
# replace immutable materials with closest mutable ones
index = ndimage.morphology.distance_transform_edt(np.in1d(material,options.immutable).reshape(grid),
return_distances = False,
return_indices = True)
material = material[index[0],
index[1],
index[2]]
immutable = np.zeros(material.shape, dtype=np.bool)
# find locations where immutable materials have been in original structure
for micro in options.immutable:
immutable += material_original == micro
# undo any changes involving immutable materials
material = np.where(immutable, material_original,material)
damask.Grid(material = material[0:grid_original[0],0:grid_original[1],0:grid_original[2]],
size = geom.size,
origin = geom.origin,
comments = geom.comments + [scriptID + ' ' + ' '.join(sys.argv[1:])],
)\
.save(sys.stdout if name is None else name)

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha5-518-g4fa97b9a3
v3.0.0-alpha5-608-g3e8d1a60d

60
python/damask/_colormap.py
View File

@ -47,27 +47,32 @@ class Colormap(mpl.colors.ListedColormap):
"""
def __eq__(self, other: object) -> bool:
def __eq__(self,
other: object) -> bool:
"""Test equality of colormaps."""
if not isinstance(other, Colormap):
return NotImplemented
return len(self.colors) == len(other.colors) \
and bool(np.all(self.colors == other.colors))
def __add__(self, other: 'Colormap') -> 'Colormap':
def __add__(self,
other: 'Colormap') -> 'Colormap':
"""Concatenate."""
return Colormap(np.vstack((self.colors,other.colors)),
f'{self.name}+{other.name}')
def __iadd__(self, other: 'Colormap') -> 'Colormap':
def __iadd__(self,
other: 'Colormap') -> 'Colormap':
"""Concatenate (in-place)."""
return self.__add__(other)
def __mul__(self, factor: int) -> 'Colormap':
def __mul__(self,
factor: int) -> 'Colormap':
"""Repeat."""
return Colormap(np.vstack([self.colors]*factor),f'{self.name}*{factor}')
def __imul__(self, factor: int) -> 'Colormap':
def __imul__(self,
factor: int) -> 'Colormap':
"""Repeat (in-place)."""
return self.__mul__(factor)
@ -161,7 +166,8 @@ class Colormap(mpl.colors.ListedColormap):
@staticmethod
def from_predefined(name: str, N: int = 256) -> 'Colormap':
def from_predefined(name: str,
N: int = 256) -> 'Colormap':
"""
Select from a set of predefined colormaps.
@ -260,10 +266,8 @@ class Colormap(mpl.colors.ListedColormap):
l,r = (field[mask].min(),field[mask].max()) if bounds is None else \
(bounds[0],bounds[1])
delta,avg = r-l,0.5*abs(r+l)
if abs(delta) * 1e8 <= avg: # delta is similar to numerical noise
l,r = l-0.5*avg*np.sign(delta),r+0.5*avg*np.sign(delta), # extend range to have actual data centered within
if abs(delta := r-l) * 1e8 <= (avg := 0.5*abs(r+l)): # delta is similar to numerical noise
l,r = (l-0.5*avg*np.sign(delta),r+0.5*avg*np.sign(delta)) # extend range to have actual data centered within
return Image.fromarray(
(np.dstack((
@ -275,7 +279,8 @@ class Colormap(mpl.colors.ListedColormap):
mode='RGBA')
def reversed(self, name: str = None) -> 'Colormap':
def reversed(self,
name: str = None) -> 'Colormap':
"""
Reverse.
@ -296,7 +301,7 @@ class Colormap(mpl.colors.ListedColormap):
>>> damask.Colormap.from_predefined('stress').reversed()
"""
rev = super(Colormap,self).reversed(name)
rev = super().reversed(name)
return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name)
@ -328,7 +333,8 @@ class Colormap(mpl.colors.ListedColormap):
return fname
def save_paraview(self, fname: FileHandle = None):
def save_paraview(self,
fname: FileHandle = None):
"""
Save as JSON file for use in Paraview.
@ -355,7 +361,8 @@ class Colormap(mpl.colors.ListedColormap):
fhandle.write('\n')
def save_ASCII(self, fname: FileHandle = None):
def save_ASCII(self,
fname: FileHandle = None):
"""
Save as ASCII file.
@ -390,7 +397,8 @@ class Colormap(mpl.colors.ListedColormap):
self._get_file_handle(fname,'.legend').write(GOM_str)
def save_gmsh(self, fname: FileHandle = None):
def save_gmsh(self,
fname: FileHandle = None):
"""
Save as ASCII file for use in gmsh.
@ -428,11 +436,11 @@ class Colormap(mpl.colors.ListedColormap):
def adjust_hue(msh_sat, msh_unsat):
"""If saturation of one of the two colors is much less than the other, hue of the less."""
if msh_sat[0] >= msh_unsat[0]:
return msh_sat[2]
else:
hSpin = msh_sat[1]/np.sin(msh_sat[1])*np.sqrt(msh_unsat[0]**2.0-msh_sat[0]**2)/msh_sat[0]
if msh_sat[2] < - np.pi/3.0: hSpin *= -1.0
return msh_sat[2] + hSpin
return msh_sat[2]
hSpin = msh_sat[1]/np.sin(msh_sat[1])*np.sqrt(msh_unsat[0]**2.0-msh_sat[0]**2)/msh_sat[0]
if msh_sat[2] < - np.pi/3.0: hSpin *= -1.0
return msh_sat[2] + hSpin
lo = np.array(low)
hi = np.array(high)
@ -445,9 +453,9 @@ class Colormap(mpl.colors.ListedColormap):
else:
lo = np.array([M_mid,0.0,0.0])
frac = 2.0*frac - 1.0
if lo[1] < 0.05 and hi[1] > 0.05:
if lo[1] < 0.05 < hi[1]:
lo[2] = adjust_hue(hi,lo)
elif lo[1] > 0.05 and hi[1] < 0.05:
elif hi[1] < 0.05 < lo[1]:
hi[2] = adjust_hue(lo,hi)
return (1.0 - frac) * lo + frac * hi
@ -476,7 +484,7 @@ class Colormap(mpl.colors.ListedColormap):
'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'brg',
'gist_rainbow', 'rainbow', 'jet', 'nipy_spectral', 'gist_ncar']}
_predefined_DAMASK = {'orientation': {'low': [0.933334,0.878432,0.878431],
_predefined_DAMASK = {'orientation': {'low': [0.933334,0.878432,0.878431], # noqa
'high': [0.250980,0.007843,0.000000]},
'strain': {'low': [0.941177,0.941177,0.870588],
'high': [0.266667,0.266667,0.000000]},
@ -621,7 +629,8 @@ class Colormap(mpl.colors.ListedColormap):
@staticmethod
def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
def _lab2xyz(lab: np.ndarray,
ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
"""
CIE Lab to CIE Xyz.
@ -652,7 +661,8 @@ class Colormap(mpl.colors.ListedColormap):
])*ref_white
@staticmethod
def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
def _xyz2lab(xyz: np.ndarray,
ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
"""
CIE Xyz to CIE Lab.

55
python/damask/_config.py
View File

@ -2,25 +2,34 @@ import copy
from io import StringIO
from collections.abc import Iterable
import abc
from pathlib import Path
from typing import Union, Dict, Any, Type, TypeVar
import numpy as np
import yaml
from ._typehints import FileHandle
from . import Rotation
MyType = TypeVar('MyType', bound='Config')
class NiceDumper(yaml.SafeDumper):
"""Make YAML readable for humans."""
def write_line_break(self, data=None):
def write_line_break(self,
data: str = None):
super().write_line_break(data)
if len(self.indents) == 1:
super().write_line_break()
def increase_indent(self, flow=False, indentless=False):
def increase_indent(self,
flow: bool = False,
indentless: bool = False):
return super().increase_indent(flow, False)
def represent_data(self, data):
def represent_data(self,
data: Any):
"""Cast Config objects and its subclasses to dict."""
if isinstance(data, dict) and type(data) != dict:
return self.represent_data(dict(data))
@ -31,14 +40,17 @@ class NiceDumper(yaml.SafeDumper):
else:
return super().represent_data(data)
def ignore_aliases(self, data):
def ignore_aliases(self,
data: Any) -> bool:
"""Do not use references to existing objects."""
return True
class Config(dict):
"""YAML-based configuration."""
def __init__(self,yml=None,**kwargs):
def __init__(self,
yml: Union[str, Dict[str, Any]] = None,
**kwargs):
"""Initialize from YAML, dict, or key=value pairs."""
if isinstance(yml,str):
kwargs.update(yaml.safe_load(yml))
@ -47,7 +59,7 @@ class Config(dict):
super().__init__(**kwargs)
def __repr__(self):
def __repr__(self) -> str:
"""Show as in file."""
output = StringIO()
self.save(output)
@ -55,14 +67,15 @@ class Config(dict):
return ''.join(output.readlines())
def __copy__(self):
def __copy__(self: MyType) -> MyType:
"""Create deep copy."""
return copy.deepcopy(self)
copy = __copy__
def __or__(self,other):
def __or__(self: MyType,
other) -> MyType:
"""
Update configuration with contents of other.
@ -76,18 +89,24 @@ class Config(dict):
updated : damask.Config
Updated configuration.
Note
----
This functionality is a backport for Python 3.8
"""
duplicate = self.copy()
duplicate.update(other)
return duplicate
def __ior__(self,other):
def __ior__(self: MyType,
other) -> MyType:
"""Update configuration with contents of other."""
return self.__or__(other)
def delete(self,keys):
def delete(self: MyType,
keys: Union[Iterable, str]) -> MyType:
"""
Remove configuration keys.
@ -109,7 +128,8 @@ class Config(dict):
@classmethod
def load(cls,fname):
def load(cls: Type[MyType],
fname: FileHandle) -> MyType:
"""
Load from yaml file.
@ -124,14 +144,15 @@ class Config(dict):
Configuration from file.
"""
try:
if isinstance(fname, (str, Path)):
fhandle = open(fname)
except TypeError:
else:
fhandle = fname
return cls(yaml.safe_load(fhandle))
def save(self,fname,**kwargs):
def save(self,
fname: FileHandle,
**kwargs):
"""
Save to yaml file.
@ -143,9 +164,9 @@ class Config(dict):
Keyword arguments parsed to yaml.dump.
"""
try:
if isinstance(fname, (str, Path)):
fhandle = open(fname,'w',newline='\n')
except TypeError:
else:
fhandle = fname
if 'width' not in kwargs:

54
python/damask/_configmaterial.py
View File

@ -1,10 +1,14 @@
import numpy as np
import h5py
from typing import Sequence, Dict, Any, Collection
from ._typehints import FileHandle
from . import Config
from . import Rotation
from . import Orientation
from . import util
from . import Table
class ConfigMaterial(Config):
"""
@ -17,7 +21,9 @@ class ConfigMaterial(Config):
"""
def __init__(self,d=None,**kwargs):
def __init__(self,
d: Dict[str, Any] = None,
**kwargs):
"""
New material configuration.
@ -30,14 +36,17 @@ class ConfigMaterial(Config):
Initial content specified as pairs of key=value.
"""
default: Collection
if d is None:
for section,default in {'material':[],'homogenization':{},'phase':{}}.items():
for section, default in {'material':[],'homogenization':{},'phase':{}}.items():
if section not in kwargs: kwargs.update({section:default})
super().__init__(d,**kwargs)
def save(self,fname='material.yaml',**kwargs):
def save(self,
fname: FileHandle = 'material.yaml',
**kwargs):
"""
Save to yaml file.
@ -53,7 +62,8 @@ class ConfigMaterial(Config):
@classmethod
def load(cls,fname='material.yaml'):
def load(cls,
fname: FileHandle = 'material.yaml') -> 'ConfigMaterial':
"""
Load from yaml file.
@ -72,10 +82,14 @@ class ConfigMaterial(Config):
@staticmethod
def load_DREAM3D(fname,
grain_data=None,cell_data=None,cell_ensemble_data='CellEnsembleData',
phases='Phases',Euler_angles='EulerAngles',phase_names='PhaseName',
base_group=None):
def load_DREAM3D(fname: str,
grain_data: str = None,
cell_data: str = None,
cell_ensemble_data: str = 'CellEnsembleData',
phases: str = 'Phases',
Euler_angles: str = 'EulerAngles',
phase_names: str = 'PhaseName',
base_group: str = None) -> 'ConfigMaterial':
"""
Load DREAM.3D (HDF5) file.
@ -154,7 +168,8 @@ class ConfigMaterial(Config):
@staticmethod
def from_table(table,**kwargs):
def from_table(table: Table,
**kwargs) -> 'ConfigMaterial':
"""
Generate from an ASCII table.
@ -207,7 +222,7 @@ class ConfigMaterial(Config):
@property
def is_complete(self):
def is_complete(self) -> bool:
"""
Check for completeness.
@ -267,12 +282,11 @@ class ConfigMaterial(Config):
if homogenization - set(self['homogenization']):
print(f'Homogenization(s) {homogenization-set(self["homogenization"])} missing')
ok = False
return ok
@property
def is_valid(self):
def is_valid(self) -> bool:
"""
Check for valid content.
@ -316,7 +330,10 @@ class ConfigMaterial(Config):
return ok
def material_rename_phase(self,mapping,ID=None,constituent=None):
def material_rename_phase(self,
mapping: Dict[str, str],
ID: Sequence[int] = None,
constituent: Sequence[int] = None) -> 'ConfigMaterial':
"""
Change phase name in material.
@ -347,7 +364,9 @@ class ConfigMaterial(Config):
return dup
def material_rename_homogenization(self,mapping,ID=None):
def material_rename_homogenization(self,
mapping: Dict[str, str],
ID: Sequence[int] = None) -> 'ConfigMaterial':
"""
Change homogenization name in material.
@ -374,7 +393,8 @@ class ConfigMaterial(Config):
return dup
def material_add(self,**kwargs):
def material_add(self,
**kwargs: Any) -> 'ConfigMaterial':
"""
Add material entries.
@ -453,7 +473,7 @@ class ConfigMaterial(Config):
N = max(N,s[0]) if len(s)>0 else N
n = max(n,s[1]) if len(s)>1 else n
mat = [{'constituents':[{} for _ in range(n)]} for _ in range(N)]
mat: Sequence[dict] = [{'constituents':[{} for _ in range(n)]} for _ in range(N)]
if 'v' not in kwargs:
shaped['v'] = np.broadcast_to(1/n,(N,n))
@ -461,7 +481,7 @@ class ConfigMaterial(Config):
map_shape = {'O':(N,n,4),'F_i':(N,n,3,3)}
for k,v in shaped.items():
target = map_shape.get(k,(N,n))
obj = np.broadcast_to(v.reshape(util.shapeshifter(v.shape,target,mode='right')),target)
obj = np.broadcast_to(v.reshape(util.shapeshifter(v.shape, target, mode = 'right')), target)
for i in range(N):
if k in ['phase','O','v','F_i']:
for j in range(n):

35
python/damask/_crystal.py
View File

@ -33,9 +33,9 @@ class Crystal():
def __init__(self,*,
family = None,
lattice = None,
a = None,b = None,c = None,
alpha = None,beta = None,gamma = None,
degrees = False):
a: float = None, b: float = None, c: float = None,
alpha: float = None, beta: float = None, gamma: float = None,
degrees: bool = False):
"""
Representation of crystal in terms of crystal family or Bravais lattice.
@ -62,7 +62,7 @@ class Crystal():
Angles are given in degrees. Defaults to False.
"""
if family not in [None] + list(lattice_symmetries.values()):
if family is not None and family not in list(lattice_symmetries.values()):
raise KeyError(f'invalid crystal family "{family}"')
if lattice is not None and family is not None and family != lattice_symmetries[lattice]:
raise KeyError(f'incompatible family "{family}" for lattice "{lattice}"')
@ -107,9 +107,6 @@ class Crystal():
if np.any([np.roll([self.alpha,self.beta,self.gamma],r)[0]
>= np.sum(np.roll([self.alpha,self.beta,self.gamma],r)[1:]) for r in range(3)]):
raise ValueError ('each lattice angle must be less than sum of others')
else:
self.a = self.b = self.c = None
self.alpha = self.beta = self.gamma = None
def __repr__(self):
@ -122,7 +119,8 @@ class Crystal():
'α={:.5g}°, β={:.5g}°, γ={:.5g}°'.format(*np.degrees(self.parameters[3:]))])
def __eq__(self,other):
def __eq__(self,
other: object) -> bool:
"""
Equal to other.
@ -132,6 +130,8 @@ class Crystal():
Crystal to check for equality.
"""
if not isinstance(other, Crystal):
return NotImplemented
return self.lattice == other.lattice and \
self.parameters == other.parameters and \
self.family == other.family
@ -139,8 +139,7 @@ class Crystal():
@property
def parameters(self):
"""Return lattice parameters a, b, c, alpha, beta, gamma."""
return (self.a,self.b,self.c,self.alpha,self.beta,self.gamma)
if hasattr(self,'a'): return (self.a,self.b,self.c,self.alpha,self.beta,self.gamma)
@property
def immutable(self):
@ -269,7 +268,7 @@ class Crystal():
https://doi.org/10.1063/1.1661333
"""
if None in self.parameters:
if self.parameters is None:
raise KeyError('missing crystal lattice parameters')
return np.array([
[1,0,0],
@ -315,7 +314,9 @@ class Crystal():
+ _lattice_points.get(self.lattice if self.lattice == 'hP' else \
self.lattice[-1],None),dtype=float)
def to_lattice(self, *, direction: np.ndarray = None, plane: np.ndarray = None) -> np.ndarray:
def to_lattice(self, *,
direction: np.ndarray = None,
plane: np.ndarray = None) -> np.ndarray:
"""
Calculate lattice vector corresponding to crystal frame direction or plane normal.
@ -339,7 +340,9 @@ class Crystal():
return np.einsum('il,...l',basis,axis)
def to_frame(self, *, uvw: np.ndarray = None, hkl: np.ndarray = None) -> np.ndarray:
def to_frame(self, *,
uvw: np.ndarray = None,
hkl: np.ndarray = None) -> np.ndarray:
"""
Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl).
@ -362,7 +365,8 @@ class Crystal():
return np.einsum('il,...l',basis,axis)
def kinematics(self, mode: str) -> Dict[str, List[np.ndarray]]:
def kinematics(self,
mode: str) -> Dict[str, List[np.ndarray]]:
"""
Return crystal kinematics systems.
@ -621,7 +625,8 @@ class Crystal():
'plane': [m[:,3:6] for m in master]}
def relation_operations(self, model: str) -> Tuple[str, Rotation]:
def relation_operations(self,
model: str) -> Tuple[str, Rotation]:
"""
Crystallographic orientation relationships for phase transformations.

104
python/damask/_grid.py
View File

@ -4,7 +4,7 @@ import warnings
import multiprocessing as mp
from functools import partial
import typing
from typing import Union, Optional, TextIO, List, Sequence, Literal
from typing import Union, Optional, TextIO, List, Sequence
from pathlib import Path
import numpy as np
@ -33,7 +33,7 @@ class Grid:
material: np.ndarray,
size: FloatSequence,
origin: FloatSequence = np.zeros(3),
comments: Union[str, Sequence[str]] = []):
comments: Union[str, Sequence[str]] = None):
"""
New geometry definition for grid solvers.
@ -53,7 +53,7 @@ class Grid:
self.material = material
self.size = size # type: ignore
self.origin = origin # type: ignore
self.comments = comments # type: ignore
self.comments = [] if comments is None else comments # type: ignore
def __repr__(self) -> str:
@ -62,8 +62,8 @@ class Grid:
mat_max = np.nanmax(self.material)
mat_N = self.N_materials
return util.srepr([
f'cells : {util.srepr(self.cells, " x ")}',
f'size : {util.srepr(self.size, " x ")} / m³',
f'cells: {util.srepr(self.cells, " × ")}',
f'size: {util.srepr(self.size, " × ")} / m³',
f'origin: {util.srepr(self.origin," ")} / m',
f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else
f' (min: {mat_min}, max: {mat_max})')
@ -77,7 +77,8 @@ class Grid:
copy = __copy__
def __eq__(self, other: object) -> bool:
def __eq__(self,
other: object) -> bool:
"""
Test equality of other.
@ -101,17 +102,18 @@ class Grid:
return self._material
@material.setter
def material(self, material: np.ndarray):
def material(self,
material: np.ndarray):
if len(material.shape) != 3:
raise ValueError(f'invalid material shape {material.shape}')
elif material.dtype not in np.sctypes['float'] and material.dtype not in np.sctypes['int']:
if material.dtype not in np.sctypes['float'] and material.dtype not in np.sctypes['int']:
raise TypeError(f'invalid material data type {material.dtype}')
else:
self._material = np.copy(material)
if self.material.dtype in np.sctypes['float'] and \
np.all(self.material == self.material.astype(int).astype(float)):
self._material = self.material.astype(int)
self._material = np.copy(material)
if self.material.dtype in np.sctypes['float'] and \
np.all(self.material == self.material.astype(int).astype(float)):
self._material = self.material.astype(int)
@property
@ -120,11 +122,12 @@ class Grid:
return self._size
@size.setter
def size(self, size: FloatSequence):
def size(self,
size: FloatSequence):
if len(size) != 3 or any(np.array(size) < 0):
raise ValueError(f'invalid size {size}')
else:
self._size = np.array(size)
self._size = np.array(size)
@property
def origin(self) -> np.ndarray:
@ -132,11 +135,12 @@ class Grid:
return self._origin
@origin.setter
def origin(self, origin: FloatSequence):
def origin(self,
origin: FloatSequence):
if len(origin) != 3:
raise ValueError(f'invalid origin {origin}')
else:
self._origin = np.array(origin)
self._origin = np.array(origin)
@property
def comments(self) -> List[str]:
@ -144,7 +148,8 @@ class Grid:
return self._comments
@comments.setter
def comments(self, comments: Union[str, Sequence[str]]):
def comments(self,
comments: Union[str, Sequence[str]]):
self._comments = [str(c) for c in comments] if isinstance(comments,list) else [str(comments)]
@ -229,8 +234,7 @@ class Grid:
content = f.readlines()
for i,line in enumerate(content[:header_length]):
items = line.split('#')[0].lower().strip().split()
key = items[0] if items else ''
if key == 'grid':
if (key := items[0] if items else '') == 'grid':
cells = np.array([ int(dict(zip(items[1::2],items[2::2]))[i]) for i in ['a','b','c']])
elif key == 'size':
size = np.array([float(dict(zip(items[1::2],items[2::2]))[i]) for i in ['x','y','z']])
@ -242,8 +246,7 @@ class Grid:
material = np.empty(int(cells.prod())) # initialize as flat array
i = 0
for line in content[header_length:]:
items = line.split('#')[0].split()
if len(items) == 3:
if len(items := line.split('#')[0].split()) == 3:
if items[1].lower() == 'of':
material_entry = np.ones(int(items[0]))*float(items[2])
elif items[1].lower() == 'to':
@ -387,7 +390,9 @@ class Grid:
@staticmethod
def _find_closest_seed(seeds: np.ndarray, weights: np.ndarray, point: np.ndarray) -> np.integer:
def _find_closest_seed(seeds: np.ndarray,
weights: np.ndarray,
point: np.ndarray) -> np.integer:
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
@staticmethod
@ -624,7 +629,9 @@ class Grid:
)
def save(self, fname: Union[str, Path], compress: bool = True):
def save(self,
fname: Union[str, Path],
compress: bool = True):
"""
Save as VTK image data file.
@ -643,7 +650,8 @@ class Grid:
v.save(fname,parallel=False,compress=compress)
def save_ASCII(self, fname: Union[str, TextIO]):
def save_ASCII(self,
fname: Union[str, TextIO]):
"""
Save as geom file.
@ -770,15 +778,16 @@ class Grid:
)
def mirror(self, directions: Sequence[str], reflect: bool = False) -> 'Grid':
def mirror(self,
directions: Sequence[str],
reflect: bool = False) -> 'Grid':
"""
Mirror grid along given directions.
Parameters
----------
directions : (sequence of) str
directions : (sequence of) {'x', 'y', 'z'}
Direction(s) along which the grid is mirrored.
Valid entries are 'x', 'y', 'z'.
reflect : bool, optional
Reflect (include) outermost layers. Defaults to False.
@ -801,8 +810,7 @@ class Grid:
# materials: 1
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):
if not set(directions).issubset(valid := ['x', 'y', 'z']):
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
limits: Sequence[Optional[int]] = [None,None] if reflect else [-2,0]
@ -822,15 +830,15 @@ class Grid:
)
def flip(self, directions: Union[Literal['x', 'y', 'z'], Sequence[Literal['x', 'y', 'z']]]) -> 'Grid':
def flip(self,
directions: Sequence[str]) -> 'Grid':
"""
Flip grid along given directions.
Parameters
----------
directions : (sequence of) str
directions : (sequence of) {'x', 'y', 'z'}
Direction(s) along which the grid is flipped.
Valid entries are 'x', 'y', 'z'.
Returns
-------
@ -838,8 +846,7 @@ class Grid:
Updated grid-based geometry.
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):
if not set(directions).issubset(valid := ['x', 'y', 'z']):
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
@ -852,7 +859,9 @@ class Grid:
)
def scale(self, cells: IntSequence, periodic: bool = True) -> 'Grid':
def scale(self,
cells: IntSequence,
periodic: bool = True) -> 'Grid':
"""
Scale grid to new cells.
@ -958,7 +967,9 @@ class Grid:
)
def rotate(self, R: Rotation, fill: int = None) -> 'Grid':
def rotate(self,
R: Rotation,
fill: int = None) -> 'Grid':
"""
Rotate grid (pad if required).
@ -1049,7 +1060,9 @@ class Grid:
)
def substitute(self, from_material: IntSequence, to_material: IntSequence) -> 'Grid':
def substitute(self,
from_material: IntSequence,
to_material: IntSequence) -> 'Grid':
"""
Substitute material indices.
@ -1150,7 +1163,9 @@ class Grid:
)
def get_grain_boundaries(self, periodic: bool = True, directions: Sequence[str] = 'xyz'):
def get_grain_boundaries(self,
periodic: bool = True,
directions: Sequence[str] = 'xyz') -> VTK:
"""
Create VTK unstructured grid containing grain boundaries.
@ -1158,9 +1173,9 @@ class Grid:
----------
periodic : bool, optional
Assume grid to be periodic. Defaults to True.
directions : (sequence of) string, optional
directions : (sequence of) {'x', 'y', 'z'}, optional
Direction(s) along which the boundaries are determined.
Valid entries are 'x', 'y', 'z'. Defaults to 'xyz'.
Defaults to 'xyz'.
Returns
-------
@ -1168,8 +1183,7 @@ class Grid:
VTK-based geometry of grain boundary network.
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):
if not set(directions).issubset(valid := ['x', 'y', 'z']):
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
o = [[0, self.cells[0]+1, np.prod(self.cells[:2]+1)+self.cells[0]+1, np.prod(self.cells[:2]+1)],

2
python/damask/_result.py
View File

@ -5,7 +5,7 @@ import os
import copy
import datetime
import warnings
import xml.etree.ElementTree as ET
import xml.etree.ElementTree as ET # noqa
import xml.dom.minidom
from pathlib import Path
from functools import partial

85
python/damask/_table.py
View File

@ -1,7 +1,7 @@
import re
import copy
from pathlib import Path
from typing import Union, Optional, Tuple, List
from typing import Union, Tuple, List
import pandas as pd
import numpy as np
@ -12,7 +12,10 @@ from . import util
class Table:
"""Manipulate multi-dimensional spreadsheet-like data."""
def __init__(self, data: np.ndarray, shapes: dict, comments: Optional[Union[str, list]] = None):
def __init__(self,
data: np.ndarray,
shapes: dict,
comments: Union[str, list] = None):
"""
New spreadsheet.
@ -41,7 +44,8 @@ class Table:
return '\n'.join(['# '+c for c in self.comments])+'\n'+data_repr
def __getitem__(self, item: Union[slice, Tuple[slice, ...]]) -> 'Table':
def __getitem__(self,
item: Union[slice, Tuple[slice, ...]]) -> 'Table':
"""
Slice the Table according to item.
@ -100,7 +104,9 @@ class Table:
copy = __copy__
def _label(self, what: Union[str, List[str]], how: str) -> List[str]:
def _label(self,
what: Union[str, List[str]],
how: str) -> List[str]:
"""
Expand labels according to data shape.
@ -131,7 +137,8 @@ class Table:
return labels
def _relabel(self, how: str):
def _relabel(self,
how: str):
"""
Modify labeling of data in-place.
@ -147,7 +154,10 @@ class Table:
self.data.columns = self._label(self.shapes,how) #type: ignore
def _add_comment(self, label: str, shape: Tuple[int, ...], info: Optional[str]):
def _add_comment(self,
label: str,
shape: Tuple[int, ...],
info: str = None):
if info is not None:
specific = f'{label}{" "+str(shape) if np.prod(shape,dtype=int) > 1 else ""}: {info}'
general = util.execution_stamp('Table')
@ -309,8 +319,7 @@ class Table:
data = np.loadtxt(content)
shapes = {'eu':3, 'pos':2, 'IQ':1, 'CI':1, 'ID':1, 'intensity':1, 'fit':1}
remainder = data.shape[1]-sum(shapes.values())
if remainder > 0: # 3.8 can do: if (remainder := data.shape[1]-sum(shapes.values())) > 0
if (remainder := data.shape[1]-sum(shapes.values())) > 0:
shapes['unknown'] = remainder
return Table(data,shapes,comments)
@ -321,7 +330,8 @@ class Table:
return list(self.shapes)
def get(self, label: str) -> np.ndarray:
def get(self,
label: str) -> np.ndarray:
"""
Get column data.
@ -341,7 +351,10 @@ class Table:
return data.astype(type(data.flatten()[0]))
def set(self, label: str, data: np.ndarray, info: str = None) -> 'Table':
def set(self,
label: str,
data: np.ndarray,
info: str = None) -> 'Table':
"""
Set column data.
@ -362,8 +375,7 @@ class Table:
"""
dup = self.copy()
dup._add_comment(label, data.shape[1:], info)
m = re.match(r'(.*)\[((\d+,)*(\d+))\]',label)
if m:
if m := re.match(r'(.*)\[((\d+,)*(\d+))\]',label):
key = m.group(1)
idx = np.ravel_multi_index(tuple(map(int,m.group(2).split(","))),
self.shapes[key])
@ -374,7 +386,10 @@ class Table:
return dup
def add(self, label: str, data: np.ndarray, info: str = None) -> 'Table':
def add(self,
label: str,
data: np.ndarray,
info: str = None) -> 'Table':
"""
Add column data.
@ -406,7 +421,8 @@ class Table:
return dup
def delete(self, label: str) -> 'Table':
def delete(self,
label: str) -> 'Table':
"""
Delete column data.
@ -427,7 +443,10 @@ class Table:
return dup
def rename(self, old: Union[str, List[str]], new: Union[str, List[str]], info: str = None) -> 'Table':
def rename(self,
old: Union[str, List[str]],
new: Union[str, List[str]],
info: str = None) -> 'Table':
"""
Rename column data.
@ -453,7 +472,9 @@ class Table:
return dup
def sort_by(self, labels: Union[str, List[str]], ascending: Union[bool, List[bool]] = True) -> 'Table':
def sort_by(self,
labels: Union[str, List[str]],
ascending: Union[bool, List[bool]] = True) -> 'Table':
"""
Sort table by values of given labels.
@ -472,8 +493,7 @@ class Table:
"""
labels_ = [labels] if isinstance(labels,str) else labels.copy()
for i,l in enumerate(labels_):
m = re.match(r'(.*)\[((\d+,)*(\d+))\]',l)
if m:
if m := re.match(r'(.*)\[((\d+,)*(\d+))\]',l):
idx = np.ravel_multi_index(tuple(map(int,m.group(2).split(','))),
self.shapes[m.group(1)])
labels_[i] = f'{1+idx}_{m.group(1)}'
@ -486,7 +506,8 @@ class Table:
return dup
def append(self, other: 'Table') -> 'Table':
def append(self,
other: 'Table') -> 'Table':
"""
Append other table vertically (similar to numpy.vstack).
@ -505,13 +526,14 @@ class Table:
"""
if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns):
raise KeyError('Labels or shapes or order do not match')
else:
dup = self.copy()
dup.data = dup.data.append(other.data,ignore_index=True)
return dup
dup = self.copy()
dup.data = dup.data.append(other.data,ignore_index=True)
return dup
def join(self, other: 'Table') -> 'Table':
def join(self,
other: 'Table') -> 'Table':
"""
Append other table horizontally (similar to numpy.hstack).
@ -530,15 +552,16 @@ class Table:
"""
if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]:
raise KeyError('Duplicated keys or row count mismatch')
else:
dup = self.copy()
dup.data = dup.data.join(other.data)
for key in other.shapes:
dup.shapes[key] = other.shapes[key]
return dup
dup = self.copy()
dup.data = dup.data.join(other.data)
for key in other.shapes:
dup.shapes[key] = other.shapes[key]
return dup
def save(self, fname: FileHandle):
def save(self,
fname: FileHandle):
"""
Save as plain text file.

3
python/damask/_typehints.py
View File

@ -9,3 +9,6 @@ import numpy as np
FloatSequence = Union[np.ndarray,Sequence[float]]
IntSequence = Union[np.ndarray,Sequence[int]]
FileHandle = Union[TextIO, str, Path]
NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.Generator]
# BitGenerator does not exists in older numpy versions
#NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.BitGenerator, np.random.Generator]

39
python/damask/_vtk.py
View File

@ -22,7 +22,8 @@ class VTK:
High-level interface to VTK.
"""
def __init__(self, vtk_data: vtk.vtkDataSet):
def __init__(self,
vtk_data: vtk.vtkDataSet):
"""
New spatial visualization.
@ -38,7 +39,9 @@ class VTK:
@staticmethod
def from_image_data(cells: IntSequence, size: FloatSequence, origin: FloatSequence = np.zeros(3)) -> 'VTK':
def from_image_data(cells: IntSequence,
size: FloatSequence,
origin: FloatSequence = np.zeros(3)) -> 'VTK':
"""
Create VTK of type vtk.vtkImageData.
@ -68,7 +71,9 @@ class VTK:
@staticmethod
def from_rectilinear_grid(grid: np.ndarray, size: FloatSequence, origin: FloatSequence = np.zeros(3)) -> 'VTK':
def from_rectilinear_grid(grid: np.ndarray,
size: FloatSequence,
origin: FloatSequence = np.zeros(3)) -> 'VTK':
"""
Create VTK of type vtk.vtkRectilinearGrid.
@ -100,7 +105,9 @@ class VTK:
@staticmethod
def from_unstructured_grid(nodes: np.ndarray, connectivity: np.ndarray, cell_type: str) -> 'VTK':
def from_unstructured_grid(nodes: np.ndarray,
connectivity: np.ndarray,
cell_type: str) -> 'VTK':
"""
Create VTK of type vtk.vtkUnstructuredGrid.
@ -195,8 +202,7 @@ class VTK:
"""
if not os.path.isfile(fname): # vtk has a strange error handling
raise FileNotFoundError(f'No such file: {fname}')
ext = Path(fname).suffix
if ext == '.vtk' or dataset_type is not None:
if (ext := Path(fname).suffix) == '.vtk' or dataset_type is not None:
reader = vtk.vtkGenericDataObjectReader()
reader.SetFileName(str(fname))
if dataset_type is None:
@ -238,7 +244,11 @@ class VTK:
def _write(writer):
"""Wrapper for parallel writing."""
writer.Write()
def save(self, fname: Union[str, Path], parallel: bool = True, compress: bool = True):
def save(self,
fname: Union[str, Path],
parallel: bool = True,
compress: bool = True):
"""
Save as VTK file.
@ -284,7 +294,9 @@ class VTK:
# Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data
# Needs support for damask.Table
def add(self, data: Union[np.ndarray, np.ma.MaskedArray], label: str = None):
def add(self,
data: Union[np.ndarray, np.ma.MaskedArray],
label: str = None):
"""
Add data to either cells or points.
@ -331,7 +343,8 @@ class VTK:
raise TypeError
def get(self, label: str) -> np.ndarray:
def get(self,
label: str) -> np.ndarray:
"""
Get either cell or point data.
@ -383,7 +396,8 @@ class VTK:
return []
def set_comments(self, comments: Union[str, List[str]]):
def set_comments(self,
comments: Union[str, List[str]]):
"""
Set comments.
@ -400,7 +414,8 @@ class VTK:
self.vtk_data.GetFieldData().AddArray(s)
def add_comments(self, comments: Union[str, List[str]]):
def add_comments(self,
comments: Union[str, List[str]]):
"""
Add comments.
@ -440,7 +455,7 @@ class VTK:
width = tk.winfo_screenwidth()
height = tk.winfo_screenheight()
tk.destroy()
except Exception as e:
except Exception:
width = 1024
height = 768

43
python/damask/grid_filters.py
View File

@ -20,7 +20,9 @@ import numpy as _np
from ._typehints import FloatSequence as _FloatSequence, IntSequence as _IntSequence
def _ks(size: _FloatSequence, cells: _IntSequence, first_order: bool = False) -> _np.ndarray:
def _ks(size: _FloatSequence,
cells: _IntSequence,
first_order: bool = False) -> _np.ndarray:
"""
Get wave numbers operator.
@ -47,7 +49,8 @@ def _ks(size: _FloatSequence, cells: _IntSequence, first_order: bool = False) ->
return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
def curl(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
def curl(size: _FloatSequence,
f: _np.ndarray) -> _np.ndarray:
u"""
Calculate curl of a vector or tensor field in Fourier space.
@ -78,7 +81,8 @@ def curl(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
return _np.fft.irfftn(curl_,axes=(0,1,2),s=f.shape[:3])
def divergence(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
def divergence(size: _FloatSequence,
f: _np.ndarray) -> _np.ndarray:
u"""
Calculate divergence of a vector or tensor field in Fourier space.
@ -105,7 +109,8 @@ def divergence(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
return _np.fft.irfftn(div_,axes=(0,1,2),s=f.shape[:3])
def gradient(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
def gradient(size: _FloatSequence,
f: _np.ndarray) -> _np.ndarray:
u"""
Calculate gradient of a scalar or vector field in Fourier space.
@ -163,7 +168,8 @@ def coordinates0_point(cells: _IntSequence,
axis = -1)
def displacement_fluct_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
def displacement_fluct_point(size: _FloatSequence,
F: _np.ndarray) -> _np.ndarray:
"""
Cell center displacement field from fluctuation part of the deformation gradient field.
@ -195,7 +201,8 @@ def displacement_fluct_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarra
return _np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3])
def displacement_avg_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
def displacement_avg_point(size: _FloatSequence,
F: _np.ndarray) -> _np.ndarray:
"""
Cell center displacement field from average part of the deformation gradient field.
@ -216,7 +223,8 @@ def displacement_avg_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_point(F.shape[:3],size))
def displacement_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
def displacement_point(size: _FloatSequence,
F: _np.ndarray) -> _np.ndarray:
"""
Cell center displacement field from deformation gradient field.
@ -236,7 +244,9 @@ def displacement_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
return displacement_avg_point(size,F) + displacement_fluct_point(size,F)
def coordinates_point(size: _FloatSequence, F: _np.ndarray, origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
def coordinates_point(size: _FloatSequence,
F: _np.ndarray,
origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
"""
Cell center positions.
@ -335,7 +345,8 @@ def coordinates0_node(cells: _IntSequence,
axis = -1)
def displacement_fluct_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
def displacement_fluct_node(size: _FloatSequence,
F: _np.ndarray) -> _np.ndarray:
"""
Nodal displacement field from fluctuation part of the deformation gradient field.
@ -355,7 +366,8 @@ def displacement_fluct_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray
return point_to_node(displacement_fluct_point(size,F))
def displacement_avg_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
def displacement_avg_node(size: _FloatSequence,
F: _np.ndarray) -> _np.ndarray:
"""
Nodal displacement field from average part of the deformation gradient field.
@ -376,7 +388,8 @@ def displacement_avg_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_node(F.shape[:3],size))
def displacement_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
def displacement_node(size: _FloatSequence,
F: _np.ndarray) -> _np.ndarray:
"""
Nodal displacement field from deformation gradient field.
@ -396,7 +409,9 @@ def displacement_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
return displacement_avg_node(size,F) + displacement_fluct_node(size,F)
def coordinates_node(size: _FloatSequence, F: _np.ndarray, origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
def coordinates_node(size: _FloatSequence,
F: _np.ndarray,
origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
"""
Nodal positions.
@ -526,7 +541,9 @@ def coordinates0_valid(coordinates0: _np.ndarray) -> bool:
return False
def regrid(size: _FloatSequence, F: _np.ndarray, cells: _IntSequence) -> _np.ndarray:
def regrid(size: _FloatSequence,
F: _np.ndarray,
cells: _IntSequence) -> _np.ndarray:
"""
Return mapping from coordinates in deformed configuration to a regular grid.

18
python/damask/mechanics.py
View File

@ -122,7 +122,9 @@ def rotation(T: _np.ndarray) -> _rotation.Rotation:
return _rotation.Rotation.from_matrix(_polar_decomposition(T,'R')[0])
def strain(F: _np.ndarray, t: str, m: float) -> _np.ndarray:
def strain(F: _np.ndarray,
t: str,
m: float) -> _np.ndarray:
"""
Calculate strain tensor (Seth–Hill family).
@ -162,7 +164,8 @@ def strain(F: _np.ndarray, t: str, m: float) -> _np.ndarray:
return eps
def stress_Cauchy(P: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
def stress_Cauchy(P: _np.ndarray,
F: _np.ndarray) -> _np.ndarray:
"""
Calculate the Cauchy stress (true stress).
@ -184,7 +187,8 @@ def stress_Cauchy(P: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
return _tensor.symmetric(_np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F))
def stress_second_Piola_Kirchhoff(P: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
def stress_second_Piola_Kirchhoff(P: _np.ndarray,
F: _np.ndarray) -> _np.ndarray:
"""
Calculate the second Piola-Kirchhoff stress.
@ -243,7 +247,8 @@ def stretch_right(T: _np.ndarray) -> _np.ndarray:
return _polar_decomposition(T,'U')[0]
def _polar_decomposition(T: _np.ndarray, requested: _Sequence[str]) -> tuple:
def _polar_decomposition(T: _np.ndarray,
requested: _Sequence[str]) -> tuple:
"""
Perform singular value decomposition.
@ -251,7 +256,7 @@ def _polar_decomposition(T: _np.ndarray, requested: _Sequence[str]) -> tuple:
----------
T : numpy.ndarray, shape (...,3,3)
Tensor of which the singular values are computed.
requested : iterable of str
requested : sequence of {'R', 'U', 'V'}
Requested outputs: ‘R’ for the rotation tensor,
‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
@ -273,7 +278,8 @@ def _polar_decomposition(T: _np.ndarray, requested: _Sequence[str]) -> tuple:
return tuple(output)
def _equivalent_Mises(T_sym: _np.ndarray, s: float) -> _np.ndarray:
def _equivalent_Mises(T_sym: _np.ndarray,
s: float) -> _np.ndarray:
"""
Base equation for Mises equivalent of a stress or strain tensor.

29
python/damask/seeds.py
View File

@ -1,3 +1,4 @@
"""Functionality for generation of seed points for Voronoi or Laguerre tessellation."""
from typing import Tuple as _Tuple
@ -5,13 +6,15 @@ from typing import Tuple as _Tuple
from scipy import spatial as _spatial
import numpy as _np
from ._typehints import FloatSequence as _FloatSequence, IntSequence as _IntSequence
from ._typehints import FloatSequence as _FloatSequence, IntSequence as _IntSequence, NumpyRngSeed as _NumpyRngSeed
from . import util as _util
from . import grid_filters as _grid_filters
def from_random(size: _FloatSequence, N_seeds: int, cells: _IntSequence = None,
rng_seed=None) -> _np.ndarray:
def from_random(size: _FloatSequence,
N_seeds: int,
cells: _IntSequence = None,
rng_seed: _NumpyRngSeed = None) -> _np.ndarray:
"""
Place seeds randomly in space.
@ -46,8 +49,12 @@ def from_random(size: _FloatSequence, N_seeds: int, cells: _IntSequence = None,
return coords
def from_Poisson_disc(size: _FloatSequence, N_seeds: int, N_candidates: int, distance: float,
periodic: bool = True, rng_seed=None) -> _np.ndarray:
def from_Poisson_disc(size: _FloatSequence,
N_seeds: int,
N_candidates: int,
distance: float,
periodic: bool = True,
rng_seed: _NumpyRngSeed = None) -> _np.ndarray:
"""
Place seeds according to a Poisson disc distribution.
@ -86,10 +93,9 @@ def from_Poisson_disc(size: _FloatSequence, N_seeds: int, N_candidates: int, dis
tree = _spatial.cKDTree(coords[:s],boxsize=size) if periodic else \
_spatial.cKDTree(coords[:s])
distances = tree.query(candidates)[0]
best = distances.argmax()
if distances[best] > distance: # require minimum separation
if distances.max() > distance: # require minimum separation
i = 0
coords[s] = candidates[best] # maximum separation to existing point cloud
coords[s] = candidates[distances.argmax()] # maximum separation to existing point cloud
s += 1
progress.update(s)
@ -99,8 +105,11 @@ def from_Poisson_disc(size: _FloatSequence, N_seeds: int, N_candidates: int, dis
return coords
def from_grid(grid, selection: _IntSequence = None, invert_selection: bool = False,
average: bool = False, periodic: bool = True) -> _Tuple[_np.ndarray, _np.ndarray]:
def from_grid(grid,
selection: _IntSequence = None,
invert_selection: bool = False,
average: bool = False,
periodic: bool = True) -> _Tuple[_np.ndarray, _np.ndarray]:
"""
Create seeds from grid description.

11
python/damask/tensor.py
View File

@ -45,7 +45,8 @@ def eigenvalues(T_sym: _np.ndarray) -> _np.ndarray:
return _np.linalg.eigvalsh(symmetric(T_sym))
def eigenvectors(T_sym: _np.ndarray, RHS: bool = False) -> _np.ndarray:
def eigenvectors(T_sym: _np.ndarray,
RHS: bool = False) -> _np.ndarray:
"""
Eigenvectors of a symmetric tensor.
@ -63,14 +64,14 @@ def eigenvectors(T_sym: _np.ndarray, RHS: bool = False) -> _np.ndarray:
associated eigenvalues.
"""
(u,v) = _np.linalg.eigh(symmetric(T_sym))
_,v = _np.linalg.eigh(symmetric(T_sym))
if RHS:
v[_np.linalg.det(v) < 0.0,:,2] *= -1.0
if RHS: v[_np.linalg.det(v) < 0.0,:,2] *= -1.0
return v
def spherical(T: _np.ndarray, tensor: bool = True) -> _np.ndarray:
def spherical(T: _np.ndarray,
tensor: bool = True) -> _np.ndarray:
"""
Calculate spherical part of a tensor.

99
python/damask/util.py
View File

@ -7,16 +7,16 @@ import subprocess
import shlex
import re
import fractions
import collections.abc as abc
from collections import abc
from functools import reduce
from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal, Optional
from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal, SupportsIndex, Sequence
from pathlib import Path
import numpy as np
import h5py
from . import version
from ._typehints import FloatSequence
from ._typehints import FloatSequence, NumpyRngSeed
# limit visibility
__all__=[
@ -54,7 +54,8 @@ _colors = {
####################################################################################################
# Functions
####################################################################################################
def srepr(msg, glue: str = '\n') -> str:
def srepr(msg,
glue: str = '\n') -> str:
r"""
Join items with glue string.
@ -76,7 +77,7 @@ def srepr(msg, glue: str = '\n') -> str:
hasattr(msg, '__iter__'))):
return glue.join(str(x) for x in msg)
else:
return msg if isinstance(msg,str) else repr(msg)
return msg if isinstance(msg,str) else repr(msg)
def emph(msg) -> str:
@ -148,7 +149,10 @@ def strikeout(msg) -> str:
return _colors['crossout']+srepr(msg)+_colors['end_color']
def run(cmd: str, wd: str = './', env: Dict[str, str] = None, timeout: int = None) -> Tuple[str, str]:
def run(cmd: str,
wd: str = './',
env: Dict[str, str] = None,
timeout: int = None) -> Tuple[str, str]:
"""
Run a command.
@ -226,15 +230,15 @@ def show_progress(iterable: Iterable,
"""
if isinstance(iterable,abc.Sequence):
if N_iter is None:
N = len(iterable)
else:
raise ValueError('N_iter given for sequence')
if N_iter is None:
N = len(iterable)
else:
raise ValueError('N_iter given for sequence')
else:
if N_iter is None:
raise ValueError('N_iter not given')
else:
N = N_iter
if N_iter is None:
raise ValueError('N_iter not given')
N = N_iter
if N <= 1:
for item in iterable:
@ -301,16 +305,20 @@ def project_equal_angle(vector: np.ndarray,
normalize : bool
Ensure unit length of input vector. Defaults to True.
keepdims : bool
Maintain three-dimensional output coordinates. Defaults to False.
Two-dimensional output uses right-handed frame spanned by
the next and next-next axis relative to the projection direction,
e.g. x-y when projecting along z and z-x when projecting along y.
Maintain three-dimensional output coordinates.
Defaults to False.
Returns
-------
coordinates : numpy.ndarray, shape (...,2 | 3)
Projected coordinates.
Notes
-----
Two-dimensional output uses right-handed frame spanned by
the next and next-next axis relative to the projection direction,
e.g. x-y when projecting along z and z-x when projecting along y.
Examples
--------
>>> import damask
@ -345,16 +353,21 @@ def project_equal_area(vector: np.ndarray,
normalize : bool
Ensure unit length of input vector. Defaults to True.
keepdims : bool
Maintain three-dimensional output coordinates. Defaults to False.
Two-dimensional output uses right-handed frame spanned by
the next and next-next axis relative to the projection direction,
e.g. x-y when projecting along z and z-x when projecting along y.
Maintain three-dimensional output coordinates.
Defaults to False.
Returns
-------
coordinates : numpy.ndarray, shape (...,2 | 3)
Projected coordinates.
Notes
-----
Two-dimensional output uses right-handed frame spanned by
the next and next-next axis relative to the projection direction,
e.g. x-y when projecting along z and z-x when projecting along y.
Examples
--------
>>> import damask
@ -373,14 +386,17 @@ def project_equal_area(vector: np.ndarray,
return np.roll(np.block([v[...,:2]/np.sqrt(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
-shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
def execution_stamp(class_name: str, function_name: str = None) -> str:
def execution_stamp(class_name: str,
function_name: str = None) -> str:
"""Timestamp the execution of a (function within a) class."""
now = datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
_function_name = '' if function_name is None else f'.{function_name}'
return f'damask.{class_name}{_function_name} v{version} ({now})'
def hybrid_IA(dist: np.ndarray, N: int, rng_seed = None) -> np.ndarray:
def hybrid_IA(dist: np.ndarray,
N: int,
rng_seed: NumpyRngSeed = None) -> np.ndarray:
"""
Hybrid integer approximation.
@ -411,7 +427,7 @@ def hybrid_IA(dist: np.ndarray, N: int, rng_seed = None) -> np.ndarray:
def shapeshifter(fro: Tuple[int, ...],
to: Tuple[int, ...],
mode: Literal['left','right'] = 'left',
keep_ones: bool = False) -> Tuple[Optional[int], ...]:
keep_ones: bool = False) -> Sequence[SupportsIndex]:
"""
Return dimensions that reshape 'fro' to become broadcastable to 'to'.
@ -447,7 +463,7 @@ def shapeshifter(fro: Tuple[int, ...],
"""
if not len(fro) and not len(to): return ()
if len(fro) == 0 and len(to) == 0: return ()
beg = dict(left ='(^.*\\b)',
right='(^.*?\\b)')
@ -455,8 +471,8 @@ def shapeshifter(fro: Tuple[int, ...],
right='(.*?\\b)')
end = dict(left ='(.*?$)',
right='(.*$)')
fro = (1,) if not len(fro) else fro
to = (1,) if not len(to) else to
fro = (1,) if len(fro) == 0 else fro
to = (1,) if len(to) == 0 else to
try:
match = re.match(beg[mode]
+f',{sep[mode]}'.join(map(lambda x: f'{x}'
@ -467,13 +483,14 @@ def shapeshifter(fro: Tuple[int, ...],
grp = match.groups()
except AssertionError:
raise ValueError(f'Shapes can not be shifted {fro} --> {to}')
fill: Tuple[Optional[int], ...] = ()
fill: Any = ()
for g,d in zip(grp,fro+(None,)):
fill += (1,)*g.count(',')+(d,)
return fill[:-1]
def shapeblender(a: Tuple[int, ...], b: Tuple[int, ...]) -> Tuple[int, ...]:
def shapeblender(a: Tuple[int, ...],
b: Tuple[int, ...]) -> Sequence[SupportsIndex]:
"""
Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
@ -517,7 +534,8 @@ def extend_docstring(extra_docstring: str) -> Callable:
return _decorator
def extended_docstring(f: Callable, extra_docstring: str) -> Callable:
def extended_docstring(f: Callable,
extra_docstring: str) -> Callable:
"""
Decorator: Combine another function's docstring with a given docstring.
@ -593,7 +611,9 @@ def DREAM3D_cell_data_group(fname: Union[str, Path]) -> str:
return cell_data_group
def Bravais_to_Miller(*, uvtw: np.ndarray = None, hkil: np.ndarray = None) -> np.ndarray:
def Bravais_to_Miller(*,
uvtw: np.ndarray = None,
hkil: np.ndarray = None) -> np.ndarray:
"""
Transform 4 Miller–Bravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl).
@ -620,7 +640,9 @@ def Bravais_to_Miller(*, uvtw: np.ndarray = None, hkil: np.ndarray = None) -> np
return np.einsum('il,...l',basis,axis)
def Miller_to_Bravais(*, uvw: np.ndarray = None, hkl: np.ndarray = None) -> np.ndarray:
def Miller_to_Bravais(*,
uvw: np.ndarray = None,
hkl: np.ndarray = None) -> np.ndarray:
"""
Transform 3 Miller indices to 4 Miller–Bravais indices of crystal direction [uvtw] or plane normal (hkil).
@ -710,7 +732,10 @@ class ProgressBar:
Works for 0-based loops, ETA is estimated by linear extrapolation.
"""
def __init__(self, total: int, prefix: str, bar_length: int):
def __init__(self,
total: int,
prefix: str,
bar_length: int):
"""
Set current time as basis for ETA estimation.
@ -733,12 +758,12 @@ class ProgressBar:
sys.stderr.write(f"{self.prefix} {'░'*self.bar_length} 0% ETA n/a")
sys.stderr.flush()
def update(self, iteration: int) -> None:
def update(self,
iteration: int) -> None:
fraction = (iteration+1) / self.total
filled_length = int(self.bar_length * fraction)
if filled_length > int(self.bar_length * self.fraction_last) or \
if filled_length := int(self.bar_length * fraction) > int(self.bar_length * self.fraction_last) or \
datetime.datetime.now() - self.time_last_update > datetime.timedelta(seconds=10):
self.time_last_update = datetime.datetime.now()
bar = '█' * filled_length + '░' * (self.bar_length - filled_length)

2
python/mypy.ini
View File

@ -1,3 +1,5 @@
[mypy]
warn_redundant_casts = True
[mypy-scipy.*]
ignore_missing_imports = True
[mypy-h5py.*]

2
python/tests/test_Orientation.py
View File

@ -287,7 +287,7 @@ class TestOrientation:
@pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('proper',[True,False])
def test_in_SST(self,family,proper):
assert Orientation(family=family).in_SST(np.zeros(3),proper)
assert Orientation(family=family).in_SST(np.zeros(3),proper)
@pytest.mark.parametrize('function',['in_SST','IPF_color'])
def test_invalid_argument(self,function):

16
python/tests/test_Result.py
View File

@ -367,13 +367,13 @@ class TestResult:
@pytest.mark.parametrize('mode',['cell','node'])
def test_coordinates(self,default,mode):
if mode == 'cell':
a = grid_filters.coordinates0_point(default.cells,default.size,default.origin)
b = default.coordinates0_point.reshape(tuple(default.cells)+(3,),order='F')
elif mode == 'node':
a = grid_filters.coordinates0_node(default.cells,default.size,default.origin)
b = default.coordinates0_node.reshape(tuple(default.cells+1)+(3,),order='F')
assert np.allclose(a,b)
if mode == 'cell':
a = grid_filters.coordinates0_point(default.cells,default.size,default.origin)
b = default.coordinates0_point.reshape(tuple(default.cells)+(3,),order='F')
elif mode == 'node':
a = grid_filters.coordinates0_node(default.cells,default.size,default.origin)
b = default.coordinates0_node.reshape(tuple(default.cells+1)+(3,),order='F')
assert np.allclose(a,b)
@pytest.mark.parametrize('output',['F','*',['P'],['P','F']],ids=range(4))
@pytest.mark.parametrize('fname',['12grains6x7x8_tensionY.hdf5'],ids=range(1))
@ -421,7 +421,7 @@ class TestResult:
def test_XDMF_datatypes(self,tmp_path,single_phase,update,ref_path):
for shape in [('scalar',()),('vector',(3,)),('tensor',(3,3)),('matrix',(12,))]:
for dtype in ['f4','f8','i1','i2','i4','i8','u1','u2','u4','u8']:
single_phase.add_calculation(f"np.ones(np.shape(#F#)[0:1]+{shape[1]},'{dtype}')",f'{shape[0]}_{dtype}')
single_phase.add_calculation(f"np.ones(np.shape(#F#)[0:1]+{shape[1]},'{dtype}')",f'{shape[0]}_{dtype}')
fname = os.path.splitext(os.path.basename(single_phase.fname))[0]+'.xdmf'
os.chdir(tmp_path)
single_phase.export_XDMF()

20
python/tests/test_Rotation.py
View File

@ -1076,19 +1076,19 @@ class TestRotation:
def test_from_fiber_component(self,N,sigma):
p = []
for run in range(5):
alpha = np.random.random()*2*np.pi,np.arccos(np.random.random())
beta = np.random.random()*2*np.pi,np.arccos(np.random.random())
alpha = np.random.random()*2*np.pi,np.arccos(np.random.random())
beta = np.random.random()*2*np.pi,np.arccos(np.random.random())
f_in_C = np.array([np.sin(alpha[0])*np.cos(alpha[1]), np.sin(alpha[0])*np.sin(alpha[1]), np.cos(alpha[0])])
f_in_S = np.array([np.sin(beta[0] )*np.cos(beta[1] ), np.sin(beta[0] )*np.sin(beta[1] ), np.cos(beta[0] )])
ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S)))
n = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0 ,normalize=True) # rotation to align fiber axis in crystal and sample system
f_in_C = np.array([np.sin(alpha[0])*np.cos(alpha[1]), np.sin(alpha[0])*np.sin(alpha[1]), np.cos(alpha[0])])
f_in_S = np.array([np.sin(beta[0] )*np.cos(beta[1] ), np.sin(beta[0] )*np.sin(beta[1] ), np.cos(beta[0] )])
ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S)))
n = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0 ,normalize=True) # rotation to align fiber axis in crystal and sample system
o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),N,False)
angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,(N,3)),n@f_in_S),-1,1))
dist = np.array(angles) * (np.random.randint(0,2,N)*2-1)
o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),N,False)
angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,(N,3)),n@f_in_S),-1,1))
dist = np.array(angles) * (np.random.randint(0,2,N)*2-1)
p.append(stats.normaltest(dist)[1])
p.append(stats.normaltest(dist)[1])
sigma_out = np.degrees(np.std(dist))
p = np.average(p)

16
python/tests/test_VTK.py
View File

@ -173,11 +173,11 @@ class TestVTK:
polyData = VTK.from_poly_data(points)
polyData.add(points,'coordinates')
if update:
polyData.save(ref_path/'polyData')
polyData.save(ref_path/'polyData')
else:
reference = VTK.load(ref_path/'polyData.vtp')
assert polyData.__repr__() == reference.__repr__() and \
np.allclose(polyData.get('coordinates'),points)
reference = VTK.load(ref_path/'polyData.vtp')
assert polyData.__repr__() == reference.__repr__() and \
np.allclose(polyData.get('coordinates'),points)
@pytest.mark.xfail(int(vtk.vtkVersion.GetVTKVersion().split('.')[0])<8, reason='missing METADATA')
def test_compare_reference_rectilinearGrid(self,update,ref_path,tmp_path):
@ -189,8 +189,8 @@ class TestVTK:
rectilinearGrid.add(np.ascontiguousarray(c),'cell')
rectilinearGrid.add(np.ascontiguousarray(n),'node')
if update:
rectilinearGrid.save(ref_path/'rectilinearGrid')
rectilinearGrid.save(ref_path/'rectilinearGrid')
else:
reference = VTK.load(ref_path/'rectilinearGrid.vtr')
assert rectilinearGrid.__repr__() == reference.__repr__() and \
np.allclose(rectilinearGrid.get('cell'),c)
reference = VTK.load(ref_path/'rectilinearGrid.vtr')
assert rectilinearGrid.__repr__() == reference.__repr__() and \
np.allclose(rectilinearGrid.get('cell'),c)

26
python/tests/test_grid_filters.py
View File

@ -8,19 +8,19 @@ from damask import seeds
class TestGridFilters:
def test_coordinates0_point(self):
size = np.random.random(3)
size = np.random.random(3) # noqa
cells = np.random.randint(8,32,(3))
coord = grid_filters.coordinates0_point(cells,size)
assert np.allclose(coord[0,0,0],size/cells*.5) and coord.shape == tuple(cells) + (3,)
def test_coordinates0_node(self):
size = np.random.random(3)
size = np.random.random(3) # noqa
cells = np.random.randint(8,32,(3))
coord = grid_filters.coordinates0_node(cells,size)
assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(cells+1) + (3,)
def test_coord0(self):
size = np.random.random(3)
size = np.random.random(3) # noqa
cells = np.random.randint(8,32,(3))
c = grid_filters.coordinates0_point(cells+1,size+size/cells)
n = grid_filters.coordinates0_node(cells,size) + size/cells*.5
@ -28,16 +28,16 @@ class TestGridFilters:
@pytest.mark.parametrize('mode',['point','node'])
def test_grid_DNA(self,mode):
"""Ensure that cellsSizeOrigin_coordinates0_xx is the inverse of coordinates0_xx."""
"""Ensure that cellsSizeOrigin_coordinates0_xx is the inverse of coordinates0_xx.""" # noqa
cells = np.random.randint(8,32,(3))
size = np.random.random(3)
origin = np.random.random(3)
coord0 = eval(f'grid_filters.coordinates0_{mode}(cells,size,origin)') # noqa
coord0 = eval(f'grid_filters.coordinates0_{mode}(cells,size,origin)') # noqa
_cells,_size,_origin = eval(f'grid_filters.cellsSizeOrigin_coordinates0_{mode}(coord0.reshape(-1,3,order="F"))')
assert np.allclose(cells,_cells) and np.allclose(size,_size) and np.allclose(origin,_origin)
def test_displacement_fluct_equivalence(self):
"""Ensure that fluctuations are periodic."""
"""Ensure that fluctuations are periodic.""" # noqa
size = np.random.random(3)
cells = np.random.randint(8,32,(3))
F = np.random.random(tuple(cells)+(3,3))
@ -45,14 +45,14 @@ class TestGridFilters:
grid_filters.point_to_node(grid_filters.displacement_fluct_point(size,F)))
def test_interpolation_to_node(self):
size = np.random.random(3)
size = np.random.random(3) # noqa
cells = np.random.randint(8,32,(3))
F = np.random.random(tuple(cells)+(3,3))
assert np.allclose(grid_filters.coordinates_node(size,F) [1:-1,1:-1,1:-1],
grid_filters.point_to_node(grid_filters.coordinates_point(size,F))[1:-1,1:-1,1:-1])
def test_interpolation_to_cell(self):
cells = np.random.randint(1,30,(3))
cells = np.random.randint(1,30,(3)) # noqa
coordinates_node_x = np.linspace(0,np.pi*2,num=cells[0]+1)
node_field_x = np.cos(coordinates_node_x)
@ -66,7 +66,7 @@ class TestGridFilters:
@pytest.mark.parametrize('mode',['point','node'])
def test_coordinates0_origin(self,mode):
origin= np.random.random(3)
origin= np.random.random(3) # noqa
size = np.random.random(3) # noqa
cells = np.random.randint(8,32,(3))
shifted = eval(f'grid_filters.coordinates0_{mode}(cells,size,origin)')
@ -79,7 +79,7 @@ class TestGridFilters:
@pytest.mark.parametrize('function',[grid_filters.displacement_avg_point,
grid_filters.displacement_avg_node])
def test_displacement_avg_vanishes(self,function):
"""Ensure that random fluctuations in F do not result in average displacement."""
"""Ensure that random fluctuations in F do not result in average displacement.""" # noqa
size = np.random.random(3)
cells = np.random.randint(8,32,(3))
F = np.random.random(tuple(cells)+(3,3))
@ -89,7 +89,7 @@ class TestGridFilters:
@pytest.mark.parametrize('function',[grid_filters.displacement_fluct_point,
grid_filters.displacement_fluct_node])
def test_displacement_fluct_vanishes(self,function):
"""Ensure that constant F does not result in fluctuating displacement."""
"""Ensure that constant F does not result in fluctuating displacement.""" # noqa
size = np.random.random(3)
cells = np.random.randint(8,32,(3))
F = np.broadcast_to(np.random.random((3,3)), tuple(cells)+(3,3))
@ -142,13 +142,13 @@ class TestGridFilters:
function(unordered,mode)
def test_regrid_identity(self):
size = np.random.random(3)
size = np.random.random(3) # noqa
cells = np.random.randint(8,32,(3))
F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
assert all(grid_filters.regrid(size,F,cells) == np.arange(cells.prod()))
def test_regrid_double_cells(self):
size = np.random.random(3)
size = np.random.random(3) # noqa
cells = np.random.randint(8,32,(3))
g = Grid.from_Voronoi_tessellation(cells,size,seeds.from_random(size,10))
F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))

3
src/CPFEM.f90
View File

@ -14,6 +14,7 @@ module CPFEM
use config
use math
use rotations
use polynomials
use lattice
use material
use phase
@ -83,12 +84,12 @@ subroutine CPFEM_initAll
call config_init
call math_init
call rotations_init
call polynomials_init
call lattice_init
call discretization_marc_init
call material_init(.false.)
call phase_init
call homogenization_init
call crystallite_init
call CPFEM_init
call config_deallocate

3
src/CPFEM2.f90
View File

@ -16,6 +16,7 @@ module CPFEM2
use config
use math
use rotations
use polynomials
use lattice
use material
use phase
@ -57,6 +58,7 @@ subroutine CPFEM_initAll
call config_init
call math_init
call rotations_init
call polynomials_init
call lattice_init
#if defined(MESH)
call discretization_mesh_init(restart=interface_restartInc>0)
@ -66,7 +68,6 @@ subroutine CPFEM_initAll
call material_init(restart=interface_restartInc>0)
call phase_init
call homogenization_init
call crystallite_init
call CPFEM_init
call config_deallocate

33
src/HDF5_utilities.f90
View File

@ -10,7 +10,8 @@ module HDF5_utilities
#include <petsc/finclude/petscsys.h>
use PETScSys
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
use MPI
use MPI_f08
use MPI, only: MPI_INFO_NULL_F90 => MPI_INFO_NULL
#endif
#endif
@ -162,9 +163,9 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel)
character, intent(in), optional :: mode
logical, intent(in), optional :: parallel
character :: m
integer(HID_T) :: plist_id
integer :: hdferr
character :: m
integer(HID_T) :: plist_id
integer :: hdferr
if (present(mode)) then
@ -178,9 +179,15 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel)
#ifdef PETSC
if (present(parallel)) then
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr)
else
call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr)
#else
if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
else
call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
#endif
end if
if(hdferr < 0) error stop 'HDF5 error'
#endif
@ -1860,7 +1867,7 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
globalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
integer, dimension(worldsize) :: &
integer(MPI_INTEGER_KIND), dimension(worldsize) :: &
readSize !< contribution of all processes
integer :: hdferr
integer(MPI_INTEGER_KIND) :: err_MPI
@ -1871,13 +1878,13 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
if(hdferr < 0) error stop 'HDF5 error'
!--------------------------------------------------------------------------------------------------
readSize = 0
readSize(worldrank+1) = int(localShape(ubound(localShape,1)))
readSize = 0_MPI_INTEGER_KIND
readSize(worldrank+1) = int(localShape(ubound(localShape,1)),MPI_INTEGER_KIND)
#ifdef PETSC
if (parallel) then
call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process
call MPI_Allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
end if
#endif
@ -1954,8 +1961,8 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
totalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id
integer, dimension(worldsize) :: writeSize !< contribution of all processes
integer(HID_T) :: dcpl
integer(MPI_INTEGER_KIND), dimension(worldsize) :: writeSize !< contribution of all processes
integer(HID_T) :: dcpl
integer :: hdferr
integer(MPI_INTEGER_KIND) :: err_MPI
integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T
@ -1974,11 +1981,11 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
!--------------------------------------------------------------------------------------------------
! determine the global data layout among all processes
writeSize = 0
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)))
writeSize = 0_MPI_INTEGER_KIND
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)),MPI_INTEGER_KIND)
#ifdef PETSC
if (parallel) then
call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process
call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
end if
#endif

4
src/IO.f90
View File

@ -584,8 +584,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
character(len=pStringLen) :: formatString
select case (warning_ID)
case (42)
msg = 'parameter has no effect'
case (47)
msg = 'no valid parameter for FFTW, using FFTW_PATIENT'
case (207)
@ -594,8 +592,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
msg = 'crystallite responds elastically'
case (601)
msg = 'stiffness close to zero'
case (700)
msg = 'unknown crystal symmetry'
case (709)
msg = 'read only the first document'
case default

6
src/YAML_parse.f90
View File

@ -191,8 +191,10 @@ logical function isScalar(line)
character(len=*), intent(in) :: line
isScalar = (.not.isKeyValue(line) .and. .not.isKey(line) .and. .not.isListItem(line) &
.and. .not.isFlow(line))
isScalar = (.not. isKeyValue(line) .and. &
.not. isKey(line) .and. &
.not. isListItem(line) .and. &
.not. isFlow(line))
end function isScalar

1
src/commercialFEM_fileList.f90
View File

@ -14,6 +14,7 @@
#include "LAPACK_interface.f90"
#include "math.f90"
#include "rotations.f90"
#include "polynomials.f90"
#include "lattice.f90"
#include "element.f90"
#include "geometry_plastic_nonlocal.f90"

6
src/constants.f90
View File

@ -9,8 +9,8 @@ module constants
public
real(pReal), parameter :: &
T_ROOM = 300.0_pReal, & !< Room temperature in K. ToDo: IUPAC: 298.15
K_B = 1.38e-23_pReal, & !< Boltzmann constant in J/Kelvin
N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol
T_ROOM = 293.15_pReal, & !< Room temperature in K (20°C)
K_B = 1.380649e-23_pReal, & !< Boltzmann constant in J/Kelvin (https://doi.org/10.1351/goldbook)
N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol (https://doi.org/10.1351/goldbook)
end module constants

11
src/grid/DAMASK_grid.f90
View File

@ -32,7 +32,7 @@ program DAMASK_grid
implicit none
type :: tLoadCase
type(rotation) :: rot !< rotation of BC
type(tRotation) :: rot !< rotation of BC
type(tBoundaryCondition) :: stress, & !< stress BC
deformation !< deformation BC (dot_F, F, or L)
real(pReal) :: t, & !< length of increment
@ -466,7 +466,14 @@ program DAMASK_grid
call MPI_Allreduce(interface_SIGUSR2,signal,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then
call mechanical_restartWrite
do field = 1, nActiveFields
select case (ID(field))
case(FIELD_MECH_ID)
call mechanical_restartWrite
case(FIELD_THERMAL_ID)
call grid_thermal_spectral_restartWrite
end select
end do
call CPFEM_restartWrite
endif
if (signal) call interface_setSIGUSR2(.false.)

144
src/grid/discretization_grid.f90
View File

@ -27,10 +27,10 @@ module discretization_grid
private
integer, dimension(3), public, protected :: &
grid !< (global) grid
cells !< (global) cells
integer, public, protected :: &
grid3, & !< (local) grid in 3rd direction
grid3Offset !< (local) grid offset in 3rd direction
cells3, & !< (local) cells in 3rd direction
cells3Offset !< (local) cells offset in 3rd direction
real(pReal), dimension(3), public, protected :: &
geomSize !< (global) physical size
real(pReal), public, protected :: &
@ -55,7 +55,7 @@ subroutine discretization_grid_init(restart)
mySize, & !< domain size of this process
origin !< (global) distance to origin
integer, dimension(3) :: &
myGrid !< domain grid of this process
myGrid !< domain cells of this process
integer, dimension(:), allocatable :: &
materialAt, materialAt_global
@ -77,7 +77,7 @@ subroutine discretization_grid_init(restart)
if (worldrank == 0) then
fileContent = IO_read(interface_geomFile)
call readVTI(grid,geomSize,origin,materialAt_global,fileContent)
call readVTI(cells,geomSize,origin,materialAt_global,fileContent)
fname = interface_geomFile
if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:)
call results_openJobFile(parallel=.false.)
@ -88,37 +88,37 @@ subroutine discretization_grid_init(restart)
end if
call MPI_Bcast(grid,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
call MPI_Bcast(cells,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1')
if (cells(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1')
call MPI_Bcast(geomSize,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Bcast(origin,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
print'(/,1x,a,3(i12,1x))', 'cells a b c: ', grid
print '(1x,a,3(es12.5,1x))', 'size x y z: ', geomSize
print '(1x,a,3(es12.5,1x))', 'origin x y z: ', origin
print'(/,1x,a,i0,a,i0,a,i0)', 'cells: ', cells(1), ' × ', cells(2), ' × ', cells(3)
print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'size: ', geomSize(1), ' × ', geomSize(2), ' × ', geomSize(3), ' / m³'
print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'origin: ', origin(1), ' ', origin(2), ' ', origin(3), ' / m'
if (worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)')
if (worldsize>cells(3)) call IO_error(894, ext_msg='number of processes exceeds cells(3)')
call fftw_mpi_init
devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), &
int(grid(2),C_INTPTR_T), &
int(grid(1),C_INTPTR_T)/2+1, &
devNull = fftw_mpi_local_size_3d(int(cells(3),C_INTPTR_T), &
int(cells(2),C_INTPTR_T), &
int(cells(1),C_INTPTR_T)/2+1, &
PETSC_COMM_WORLD, &
z, & ! domain grid size along z
z_offset) ! domain grid offset along z
z, & ! domain cells size along z
z_offset) ! domain cells offset along z
if (z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes')
grid3 = int(z)
grid3Offset = int(z_offset)
size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal)
size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal)
myGrid = [grid(1:2),grid3]
cells3 = int(z)
cells3Offset = int(z_offset)
size3 = geomSize(3)*real(cells3,pReal) /real(cells(3),pReal)
size3Offset = geomSize(3)*real(cells3Offset,pReal)/real(cells(3),pReal)
myGrid = [cells(1:2),cells3]
mySize = [geomSize(1:2),size3]
call MPI_Gather(product(grid(1:2))*grid3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,&
call MPI_Gather(product(cells(1:2))*cells3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,&
1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Gather(product(myGrid), 1_MPI_INTEGER_KIND,MPI_INTEGER,sendcounts,&
@ -131,10 +131,10 @@ subroutine discretization_grid_init(restart)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call discretization_init(materialAt, &
IPcoordinates0(myGrid,mySize,grid3Offset), &
Nodes0(myGrid,mySize,grid3Offset),&
merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write top layer...
(grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process
IPcoordinates0(myGrid,mySize,cells3Offset), &
Nodes0(myGrid,mySize,cells3Offset),&
merge((cells(1)+1) * (cells(2)+1) * (cells3+1),& ! write top layer...
(cells(1)+1) * (cells(2)+1) * cells3,& ! ...unless not last process
worldrank+1==worldsize))
!--------------------------------------------------------------------------------------------------
@ -142,7 +142,7 @@ subroutine discretization_grid_init(restart)
if (.not. restart) then
call results_openJobFile
call results_closeGroup(results_addGroup('geometry'))
call results_addAttribute('cells', grid, '/geometry')
call results_addAttribute('cells', cells, '/geometry')
call results_addAttribute('size', geomSize,'/geometry')
call results_addAttribute('origin',origin, '/geometry')
call results_closeJobFile
@ -170,13 +170,13 @@ end subroutine discretization_grid_init
!> @brief Parse vtk image data (.vti)
!> @details https://vtk.org/Wiki/VTK_XML_Formats
!--------------------------------------------------------------------------------------------------
subroutine readVTI(grid,geomSize,origin,material, &
subroutine readVTI(cells,geomSize,origin,material, &
fileContent)
integer, dimension(3), intent(out) :: &
grid ! grid (across all processes!)
cells ! cells (across all processes!)
real(pReal), dimension(3), intent(out) :: &
geomSize, & ! size (across all processes!)
geomSize, & ! size (across all processes!)
origin ! origin (across all processes!)
integer, dimension(:), intent(out), allocatable :: &
material
@ -190,7 +190,7 @@ subroutine readVTI(grid,geomSize,origin,material, &
s
grid = -1
cells = -1
geomSize = -1.0_pReal
inFile = .false.
@ -215,7 +215,7 @@ subroutine readVTI(grid,geomSize,origin,material, &
if (.not. inImage) then
if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then
inImage = .true.
call cellsSizeOrigin(grid,geomSize,origin,fileContent(startPos:endPos))
call cellsSizeOrigin(cells,geomSize,origin,fileContent(startPos:endPos))
end if
else
if (index(fileContent(startPos:endPos),'<CellData',kind=pI64) /= 0_pI64) then
@ -246,10 +246,10 @@ subroutine readVTI(grid,geomSize,origin,material, &
end do
if (.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found')
if (size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)')
if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if (any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid')
if (.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found')
if (size(material) /= product(cells)) call IO_error(error_ID = 844, ext_msg='size(material)')
if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if (any(cells<1)) call IO_error(error_ID = 844, ext_msg='cells')
material = material + 1
if (any(material<1)) call IO_error(error_ID = 844, ext_msg='material ID < 0')
@ -502,13 +502,13 @@ end subroutine readVTI
!---------------------------------------------------------------------------------------------------
!> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element)
!---------------------------------------------------------------------------------------------------
function IPcoordinates0(grid,geomSize,grid3Offset)
function IPcoordinates0(cells,geomSize,cells3Offset)
integer, dimension(3), intent(in) :: grid ! grid (for this process!)
integer, dimension(3), intent(in) :: cells ! cells (for this process!)
real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!)
integer, intent(in) :: grid3Offset ! grid(3) offset
integer, intent(in) :: cells3Offset ! cells(3) offset
real(pReal), dimension(3,product(grid)) :: ipCoordinates0
real(pReal), dimension(3,product(cells)) :: ipCoordinates0
integer :: &
a,b,c, &
@ -516,9 +516,9 @@ function IPcoordinates0(grid,geomSize,grid3Offset)
i = 0
do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1)
do c = 1, cells(3); do b = 1, cells(2); do a = 1, cells(1)
i = i + 1
IPcoordinates0(1:3,i) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal)
IPcoordinates0(1:3,i) = geomSize/real(cells,pReal) * (real([a,b,cells3Offset+c],pReal) -0.5_pReal)
end do; end do; end do
end function IPcoordinates0
@ -527,22 +527,22 @@ end function IPcoordinates0
!---------------------------------------------------------------------------------------------------
!> @brief Calculate position of undeformed nodes (pretend to be an element)
!---------------------------------------------------------------------------------------------------
pure function nodes0(grid,geomSize,grid3Offset)
pure function nodes0(cells,geomSize,cells3Offset)
integer, dimension(3), intent(in) :: grid ! grid (for this process!)
integer, dimension(3), intent(in) :: cells ! cells (for this process!)
real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!)
integer, intent(in) :: grid3Offset ! grid(3) offset
integer, intent(in) :: cells3Offset ! cells(3) offset
real(pReal), dimension(3,product(grid+1)) :: nodes0
real(pReal), dimension(3,product(cells+1)) :: nodes0
integer :: &
a,b,c, &
n
n = 0
do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1)
do c = 0, cells3; do b = 0, cells(2); do a = 0, cells(1)
n = n + 1
nodes0(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal)
nodes0(1:3,n) = geomSize/real(cells,pReal) * real([a,b,cells3Offset+c],pReal)
end do; end do; end do
end function nodes0
@ -551,17 +551,17 @@ end function nodes0
!--------------------------------------------------------------------------------------------------
!> @brief Calculate IP interface areas
!--------------------------------------------------------------------------------------------------
pure function cellSurfaceArea(geomSize,grid)
pure function cellSurfaceArea(geomSize,cells)
real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!)
integer, dimension(3), intent(in) :: grid ! grid (for this process!)
integer, dimension(3), intent(in) :: cells ! cells (for this process!)
real(pReal), dimension(6,1,product(grid)) :: cellSurfaceArea
real(pReal), dimension(6,1,product(cells)) :: cellSurfaceArea
cellSurfaceArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3))
cellSurfaceArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1))
cellSurfaceArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2))
cellSurfaceArea(1:2,1,:) = geomSize(2)/real(cells(2)) * geomSize(3)/real(cells(3))
cellSurfaceArea(3:4,1,:) = geomSize(3)/real(cells(3)) * geomSize(1)/real(cells(1))
cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1)) * geomSize(2)/real(cells(2))
end function cellSurfaceArea
@ -588,42 +588,42 @@ end function cellSurfaceNormal
!--------------------------------------------------------------------------------------------------
!> @brief Build IP neighborhood relations
!--------------------------------------------------------------------------------------------------
pure function IPneighborhood(grid)
pure function IPneighborhood(cells)
integer, dimension(3), intent(in) :: grid ! grid (for this process!)
integer, dimension(3), intent(in) :: cells ! cells (for this process!)
integer, dimension(3,6,1,product(grid)) :: IPneighborhood !< 6 neighboring IPs as [element ID, IP ID, face ID]
integer, dimension(3,6,1,product(cells)) :: IPneighborhood !< 6 neighboring IPs as [element ID, IP ID, face ID]
integer :: &
x,y,z, &
e
e = 0
do z = 0,grid(3)-1; do y = 0,grid(2)-1; do x = 0,grid(1)-1
do z = 0,cells(3)-1; do y = 0,cells(2)-1; do x = 0,cells(1)-1
e = e + 1
! element ID
IPneighborhood(1,1,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x+1,grid(1)) &
IPneighborhood(1,1,1,e) = z * cells(1) * cells(2) &
+ y * cells(1) &
+ modulo(x+1,cells(1)) &
+ 1
IPneighborhood(1,2,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x-1,grid(1)) &
IPneighborhood(1,2,1,e) = z * cells(1) * cells(2) &
+ y * cells(1) &
+ modulo(x-1,cells(1)) &
+ 1
IPneighborhood(1,3,1,e) = z * grid(1) * grid(2) &
+ modulo(y+1,grid(2)) * grid(1) &
IPneighborhood(1,3,1,e) = z * cells(1) * cells(2) &
+ modulo(y+1,cells(2)) * cells(1) &
+ x &
+ 1
IPneighborhood(1,4,1,e) = z * grid(1) * grid(2) &
+ modulo(y-1,grid(2)) * grid(1) &
IPneighborhood(1,4,1,e) = z * cells(1) * cells(2) &
+ modulo(y-1,cells(2)) * cells(1) &
+ x &
+ 1
IPneighborhood(1,5,1,e) = modulo(z+1,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
IPneighborhood(1,5,1,e) = modulo(z+1,cells(3)) * cells(1) * cells(2) &
+ y * cells(1) &
+ x &
+ 1
IPneighborhood(1,6,1,e) = modulo(z-1,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
IPneighborhood(1,6,1,e) = modulo(z-1,cells(3)) * cells(1) * cells(2) &
+ y * cells(1) &
+ x &
+ 1
! IP ID

36
src/grid/grid_damage_spectral.f90
View File

@ -106,9 +106,9 @@ subroutine grid_damage_spectral_init()
!--------------------------------------------------------------------------------------------------
! init fields
allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_current(cells(1),cells(2),cells3), source=1.0_pReal)
allocate(phi_lastInc(cells(1),cells(2),cells3), source=1.0_pReal)
allocate(phi_stagInc(cells(1),cells(2),cells3), source=1.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -117,23 +117,23 @@ subroutine grid_damage_spectral_init()
call SNESSetOptionsPrefix(SNES_damage,'damage_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3D(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
1_pPetscInt, 0_pPetscInt, & ! #dof (phi, scalar), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
damage_grid,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call DMsetFromOptions(damage_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(damage_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (cells x 1, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
@ -213,7 +213,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
!--------------------------------------------------------------------------------------------------
! updating damage state
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
call homogenization_set_phi(phi_current(i,j,k),ce)
end do; end do; end do
@ -255,7 +255,7 @@ subroutine grid_damage_spectral_forward(cutBack)
call DMDAVecRestoreArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
call homogenization_set_phi(phi_current(i,j,k),ce)
end do; end do; end do
@ -289,12 +289,12 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current
scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_current
call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k))
end do; end do; end do
@ -302,7 +302,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
+ homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
@ -315,14 +315,14 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t)
call utilities_FFTscalarBackward
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness
where(scalarField_real(1:cells(1),1:cells(2),1:cells3) > phi_lastInc) &
scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_lastInc
where(scalarField_real(1:cells(1),1:cells(2),1:cells3) < num%residualStiffness) &
scalarField_real(1:cells(1),1:cells(2),1:cells3) = num%residualStiffness
!--------------------------------------------------------------------------------------------------
! constructing residual
r = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current
r = scalarField_real(1:cells(1),1:cells(2),1:cells3) - phi_current
err_PETSc = 0
end subroutine formResidual
@ -339,7 +339,7 @@ subroutine updateReference()
K_ref = 0.0_pReal
mu_ref = 0.0_pReal
do ce = 1, product(grid(1:2))*grid3
do ce = 1, product(cells(1:2))*cells3
K_ref = K_ref + homogenization_K_phi(ce)
mu_ref = mu_ref + homogenization_mu_phi(ce)
end do

56
src/grid/grid_mech_FEM.f90
View File

@ -153,9 +153,9 @@ subroutine grid_mechanical_FEM_init
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate(F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(P_current (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(F_lastInc (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -164,16 +164,16 @@ subroutine grid_mechanical_FEM_init
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if(err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, &
DMDA_STENCIL_BOX, &
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
3_pPetscInt, 1_pPetscInt, & ! #dof (u, vector), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetFromOptions(mechanical_grid,err_PETSc)
@ -214,7 +214,7 @@ subroutine grid_mechanical_FEM_init
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,err_PETSc)
CHKERRQ(err_PETSc)
delta = geomSize/real(grid,pReal) ! grid spacing
delta = geomSize/real(cells,pReal) ! grid spacing
detJ = product(delta) ! cell volume
BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), &
@ -255,11 +255,11 @@ subroutine grid_mechanical_FEM_init
call HDF5_read(u_lastInc,groupHandle,'u_lastInc')
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity
F = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3)
endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(F)
call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F
@ -339,7 +339,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(rotation), intent(in) :: &
type(tRotation), intent(in) :: &
rotation_BC
PetscErrorCode :: err_PETSc
PetscScalar, pointer, dimension(:,:,:,:) :: &
@ -386,7 +386,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
F_lastInc = F
homogenization_F0 = reshape(F, [3,3,product(grid(1:2))*grid3])
homogenization_F0 = reshape(F, [3,3,product(cells(1:2))*cells3])
endif
!--------------------------------------------------------------------------------------------------
@ -556,13 +556,13 @@ subroutine formResidual(da_local,x_local, &
! get deformation gradient
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
ctr = 0
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
enddo; enddo; enddo
F(1:3,1:3,i,j,k-grid3offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
F(1:3,1:3,i,j,k-cells3Offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
@ -589,14 +589,14 @@ subroutine formResidual(da_local,x_local, &
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
ele = 0
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
ctr = 0
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
enddo; enddo; enddo
ele = ele + 1
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-grid3offset)))*detJ + &
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-cells3Offset)))*detJ + &
matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
@ -615,17 +615,17 @@ subroutine formResidual(da_local,x_local, &
! applying boundary conditions
call DMDAVecGetArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
if (grid3offset == 0) then
r(0:2,0, 0, 0) = 0.0_pReal
r(0:2,grid(1),0, 0) = 0.0_pReal
r(0:2,0, grid(2),0) = 0.0_pReal
r(0:2,grid(1),grid(2),0) = 0.0_pReal
if (cells3Offset == 0) then
r(0:2,0, 0, 0) = 0.0_pReal
r(0:2,cells(1),0, 0) = 0.0_pReal
r(0:2,0, cells(2),0) = 0.0_pReal
r(0:2,cells(1),cells(2),0) = 0.0_pReal
end if
if (grid3+grid3offset == grid(3)) then
r(0:2,0, 0, grid(3)) = 0.0_pReal
r(0:2,grid(1),0, grid(3)) = 0.0_pReal
r(0:2,0, grid(2),grid(3)) = 0.0_pReal
r(0:2,grid(1),grid(2),grid(3)) = 0.0_pReal
if (cells3+cells3Offset == cells(3)) then
r(0:2,0, 0, cells(3)) = 0.0_pReal
r(0:2,cells(1),0, cells(3)) = 0.0_pReal
r(0:2,0, cells(2),cells(3)) = 0.0_pReal
r(0:2,cells(1),cells(2),cells(3)) = 0.0_pReal
end if
call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
@ -663,7 +663,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
call MatZeroEntries(Jac,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
ctr = 0
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1
@ -719,7 +719,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
ce = ce + 1
x_scal(0:2,i-1,j-1,k-1) = discretization_IPcoords(1:3,ce)
enddo; enddo; enddo

42
src/grid/grid_mech_spectral_basic.f90
View File

@ -79,6 +79,12 @@ module grid_mechanical_spectral_basic
err_BC, & !< deviation from stress BC
err_div !< RMS of div of P
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
type(MPI_Status) :: status
#else
integer, dimension(MPI_STATUS_SIZE) :: status
#endif
integer :: &
totalIter = 0 !< total iteration in current increment
@ -96,7 +102,7 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_basic_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: P
PetscErrorCode :: err_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI
PetscScalar, pointer, dimension(:,:,:,:) :: &
@ -153,8 +159,8 @@ subroutine grid_mechanical_spectral_basic_init
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate(F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_lastInc(3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(Fdot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -163,23 +169,23 @@ subroutine grid_mechanical_spectral_basic_init
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
9_pPetscInt, 0_pPetscInt, & ! #dof (F, tensor), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
da,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call DMsetFromOptions(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (cells x 9, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
@ -217,11 +223,11 @@ subroutine grid_mechanical_spectral_basic_init
call HDF5_read(F_lastInc,groupHandle,'F_lastInc')
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity
F = reshape(F_lastInc,[9,cells(1),cells(2),cells3])
end if restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
@ -244,7 +250,7 @@ subroutine grid_mechanical_spectral_basic_init
call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', &
MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI)
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_close(fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@ -310,7 +316,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(rotation), intent(in) :: &
type(tRotation), intent(in) :: &
rotation_BC
PetscErrorCode :: err_PETSc
PetscScalar, pointer, dimension(:,:,:,:) :: F
@ -343,11 +349,11 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
end if
Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, &
rotation_BC%rotate(F_aimDot,active=.true.))
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
F_lastInc = reshape(F,[3,3,cells(1),cells(2),cells3])
homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3])
end if
!--------------------------------------------------------------------------------------------------
@ -359,7 +365,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
+ merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t 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,cells(1),cells(2),cells3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
@ -530,7 +536,7 @@ subroutine formResidual(in, F, &
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r ! store fPK field for subsequent FFT forward transform
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
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
@ -538,7 +544,7 @@ subroutine formResidual(in, F, &
!--------------------------------------------------------------------------------------------------
! constructing residual
r = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
r = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
end subroutine formResidual

60
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -90,6 +90,12 @@ module grid_mechanical_spectral_polarisation
err_curl, & !< RMS of curl of F
err_div !< RMS of div of P
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
type(MPI_Status) :: status
#else
integer, dimension(MPI_STATUS_SIZE) :: status
#endif
integer :: &
totalIter = 0 !< total iteration in current increment
@ -107,7 +113,7 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_polarisation_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: P
PetscErrorCode :: err_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI
PetscScalar, pointer, dimension(:,:,:,:) :: &
@ -171,10 +177,10 @@ subroutine grid_mechanical_spectral_polarisation_init
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_lastInc (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(Fdot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(F_tau_lastInc(3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
allocate(F_tauDot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -183,23 +189,23 @@ subroutine grid_mechanical_spectral_polarisation_init
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
18_pPetscInt, 0_pPetscInt, & ! #dof (2xtensor), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
da,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call DMsetFromOptions(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (cells x 18, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
@ -241,13 +247,13 @@ subroutine grid_mechanical_spectral_polarisation_init
call HDF5_read(F_tau_lastInc,groupHandle,'F_tau_lastInc')
elseif (interface_restartInc == 0) then restartRead
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity
F = reshape(F_lastInc,[9,cells(1),cells(2),cells3])
F_tau = 2.0_pReal*F
F_tau_lastInc = 2.0_pReal*F_lastInc
end if restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
@ -270,7 +276,7 @@ subroutine grid_mechanical_spectral_polarisation_init
call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', &
MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI)
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_close(fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@ -342,7 +348,7 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(rotation), intent(in) :: &
type(tRotation), intent(in) :: &
rotation_BC
PetscErrorCode :: err_PETSc
PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau
@ -379,15 +385,15 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
end if
Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, &
rotation_BC%rotate(F_aimDot,active=.true.))
F_tauDot = utilities_calculateRate(guess, &
F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), Delta_t_old, &
F_tau_lastInc,reshape(F_tau,[3,3,cells(1),cells(2),cells3]), Delta_t_old, &
rotation_BC%rotate(F_aimDot,active=.true.))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
F_lastInc = reshape(F, [3,3,cells(1),cells(2),cells3])
F_tau_lastInc = reshape(F_tau,[3,3,cells(1),cells(2),cells3])
homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3])
end if
!--------------------------------------------------------------------------------------------------
@ -400,12 +406,12 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),&
[9,grid(1),grid(2),grid3])
[9,cells(1),cells(2),cells3])
if (guess) then
F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), &
[9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition
[9,cells(1),cells(2),cells3]) ! does not have any average value as boundary condition
else
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
F_lambda33 = math_I3 &
+ math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, &
@ -597,7 +603,7 @@ subroutine formResidual(in, FandF_tau, &
!--------------------------------------------------------------------------------------------------
!
tensorField_real = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
tensorField_real(1:3,1:3,i,j,k) = &
num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
num%alpha*matmul(F(1:3,1:3,i,j,k), &
@ -612,7 +618,7 @@ subroutine formResidual(in, FandF_tau, &
!--------------------------------------------------------------------------------------------------
! constructing residual
r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
@ -629,14 +635,14 @@ subroutine formResidual(in, FandF_tau, &
params%stress_mask)))
! calculate divergence
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r_F !< stress field in disguise
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r_F !< stress field in disguise
call utilities_FFTtensorForward
err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
e = e + 1
r_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), &
@ -648,7 +654,7 @@ subroutine formResidual(in, FandF_tau, &
!--------------------------------------------------------------------------------------------------
! calculating curl
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:cells(1),1:cells(2),1:cells3) = F
call utilities_FFTtensorForward
err_curl = utilities_curlRMS()

89
src/grid/grid_thermal_spectral.f90
View File

@ -16,6 +16,9 @@ module grid_thermal_spectral
use prec
use parallelization
use IO
use DAMASK_interface
use HDF5_utilities
use HDF5
use spectral_utilities
use discretization_grid
use homogenization
@ -54,13 +57,13 @@ module grid_thermal_spectral
public :: &
grid_thermal_spectral_init, &
grid_thermal_spectral_solution, &
grid_thermal_spectral_restartWrite, &
grid_thermal_spectral_forward
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data
! ToDo: Restart not implemented
!--------------------------------------------------------------------------------------------------
subroutine grid_thermal_spectral_init(T_0)
@ -72,6 +75,7 @@ subroutine grid_thermal_spectral_init(T_0)
PetscScalar, dimension(:,:,:), pointer :: T_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI
PetscErrorCode :: err_PETSc
integer(HID_T) :: fileHandle, groupHandle
class(tNode), pointer :: &
num_grid
@ -101,15 +105,9 @@ subroutine grid_thermal_spectral_init(T_0)
!--------------------------------------------------------------------------------------------------
! init fields
allocate(T_current(grid(1),grid(2),grid3), source=T_0)
allocate(T_lastInc(grid(1),grid(2),grid3), source=T_0)
allocate(T_stagInc(grid(1),grid(2),grid3), source=T_0)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1
call homogenization_thermal_setField(T_0,0.0_pReal,ce)
end do; end do; end do
allocate(T_current(cells(1),cells(2),cells3), source=T_0)
allocate(T_lastInc(cells(1),cells(2),cells3), source=T_0)
allocate(T_stagInc(cells(1),cells(2),cells3), source=T_0)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -118,23 +116,23 @@ subroutine grid_thermal_spectral_init(T_0)
call SNESSetOptionsPrefix(SNES_thermal,'thermal_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
localK(worldrank) = int(cells3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call DMDACreate3D(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid
int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells
1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), &
1_pPetscInt, 0_pPetscInt, & ! #dof (T, scalar), ghost boundary width (domain overlap)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells
thermal_grid,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call DMsetFromOptions(thermal_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(thermal_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (cells x 1, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
@ -142,6 +140,24 @@ subroutine grid_thermal_spectral_init(T_0)
CHKERRQ(err_PETSc)
call SNESSetFromOptions(SNES_thermal,err_PETSc) ! pull it all together with additional CLI arguments
CHKERRQ(err_PETSc)
restartRead: if (interface_restartInc > 0) then
print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_read(T_current,groupHandle,'T',.false.)
call HDF5_read(T_lastInc,groupHandle,'T_lastInc',.false.)
end if restartRead
ce = 0
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)
ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),0.0_pReal,ce)
end do; end do; end do
call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
T_PETSc = T_current
@ -198,7 +214,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
!--------------------------------------------------------------------------------------------------
! updating thermal state
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
end do; end do; end do
@ -241,7 +257,7 @@ subroutine grid_thermal_spectral_forward(cutBack)
call DMDAVecRestoreArrayF90(dm_local,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
end do; end do; end do
@ -253,6 +269,37 @@ subroutine grid_thermal_spectral_forward(cutBack)
end subroutine grid_thermal_spectral_forward
!--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file
!--------------------------------------------------------------------------------------------------
subroutine grid_thermal_spectral_restartWrite
PetscErrorCode :: err_PETSc
DM :: dm_local
integer(HID_T) :: fileHandle, groupHandle
PetscScalar, dimension(:,:,:), pointer :: T
call SNESGetDM(SNES_thermal,dm_local,err_PETSc);
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(dm_local,solution_vec,T,err_PETSc);
CHKERRQ(err_PETSc)
print'(1x,a)', 'writing thermal solver data required for restart to file'; flush(IO_STDOUT)
fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','a')
groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_write(T,groupHandle,'T')
call HDF5_write(T_lastInc,groupHandle,'T_lastInc')
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
call DMDAVecRestoreArrayF90(dm_local,solution_vec,T,err_PETSc);
CHKERRQ(err_PETSc)
end subroutine grid_thermal_spectral_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief forms the spectral thermal residual vector
!--------------------------------------------------------------------------------------------------
@ -274,12 +321,12 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = T_current
scalarField_real(1:cells(1),1:cells(2),1:cells3) = T_current
call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of temperature field
call utilities_FFTvectorBackward
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k))
end do; end do; end do
@ -287,7 +334,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field
call utilities_FFTscalarBackward
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) &
+ homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
@ -302,7 +349,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
!--------------------------------------------------------------------------------------------------
! constructing residual
r = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3)
r = T_current - scalarField_real(1:cells(1),1:cells(2),1:cells3)
err_PETSc = 0
end subroutine formResidual
@ -319,7 +366,7 @@ subroutine updateReference()
K_ref = 0.0_pReal
mu_ref = 0.0_pReal
do ce = 1, product(grid(1:2))*grid3
do ce = 1, product(cells(1:2))*cells3
K_ref = K_ref + homogenization_K_T(ce)
mu_ref = mu_ref + homogenization_mu_T(ce)
end do

208
src/grid/spectral_utilities.f90
View File

@ -29,9 +29,9 @@ module spectral_utilities
include 'fftw3-mpi.f03'
!--------------------------------------------------------------------------------------------------
! grid related information information
! grid related information
real(pReal), protected, public :: wgt !< weighting factor 1/Nelems
integer, protected, public :: grid1Red !< grid(1)/2
integer, protected, public :: grid1Red !< cells(1)/2
real(pReal), protected, public, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence
!--------------------------------------------------------------------------------------------------
@ -86,7 +86,7 @@ module spectral_utilities
type, public :: tSolutionParams
real(pReal), dimension(3,3) :: stress_BC
logical, dimension(3,3) :: stress_mask
type(rotation) :: rotation_BC
type(tRotation) :: rotation_BC
real(pReal) :: Delta_t
end type tSolutionParams
@ -201,8 +201,8 @@ subroutine spectral_utilities_init
num_grid%get_asString('PETSc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc)
grid1Red = grid(1)/2 + 1
wgt = 1.0/real(product(grid),pReal)
grid1Red = cells(1)/2 + 1
wgt = 1.0/real(product(cells),pReal)
num%memory_efficient = num_grid%get_asInt('memory_efficient', defaultVal=1) > 0 ! ToDo: should be logical in YAML file
num%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2)
@ -231,9 +231,9 @@ subroutine spectral_utilities_init
enddo
elseif (num%divergence_correction == 2) then
do j = 1, 3
if ( j /= int(minloc(geomSize/real(grid,pReal),1)) &
.and. j /= int(maxloc(geomSize/real(grid,pReal),1))) &
scaledGeomSize = geomSize/geomSize(j)*real(grid(j),pReal)
if ( j /= int(minloc(geomSize/real(cells,pReal),1)) &
.and. j /= int(maxloc(geomSize/real(cells,pReal),1))) &
scaledGeomSize = geomSize/geomSize(j)*real(cells(j),pReal)
enddo
else
scaledGeomSize = geomSize
@ -262,11 +262,11 @@ subroutine spectral_utilities_init
!--------------------------------------------------------------------------------------------------
! MPI allocation
gridFFTW = int(grid,C_INTPTR_T)
gridFFTW = int(cells,C_INTPTR_T)
alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, &
PETSC_COMM_WORLD, local_K, local_K_offset)
allocate (xi1st (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension
allocate (xi2nd (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension
allocate (xi1st (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension
allocate (xi2nd (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension
tensorField = fftw_alloc_complex(tensorSize*alloc_local)
call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, &
@ -327,27 +327,27 @@ subroutine spectral_utilities_init
!--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
do k = grid3Offset+1, grid3Offset+grid3
do k = cells3Offset+1, cells3Offset+cells3
k_s(3) = k - 1
if (k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do j = 1, grid(2)
if (k > cells(3)/2 + 1) k_s(3) = k_s(3) - cells(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do j = 1, cells(2)
k_s(2) = j - 1
if (j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
if (j > cells(2)/2 + 1) k_s(2) = k_s(2) - cells(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do i = 1, grid1Red
k_s(1) = i - 1 ! symmetry, junst running from 0,1,...,N/2,N/2+1
xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s)
where(mod(grid,2)==0 .and. [i,j,k] == grid/2+1 .and. &
xi2nd(1:3,i,j,k-cells3Offset) = utilities_getFreqDerivative(k_s)
where(mod(cells,2)==0 .and. [i,j,k] == cells/2+1 .and. &
spectral_derivative_ID == DERIVATIVE_CONTINUOUS_ID) ! for even grids, set the Nyquist Freq component to 0.0
xi1st(1:3,i,j,k-grid3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal)
xi1st(1:3,i,j,k-cells3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal)
elsewhere
xi1st(1:3,i,j,k-grid3Offset) = xi2nd(1:3,i,j,k-grid3Offset)
xi1st(1:3,i,j,k-cells3Offset) = xi2nd(1:3,i,j,k-cells3Offset)
endwhere
enddo; enddo; enddo
if (num%memory_efficient) then ! allocate just single fourth order tensor
if (num%memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal))
else ! precalculation of gamma_hat field
allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
allocate (gamma_hat(3,3,3,3,grid1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
endif
end subroutine spectral_utilities_init
@ -373,10 +373,10 @@ subroutine utilities_updateGamma(C)
if (.not. num%memory_efficient) then
gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A
do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, grid1Red
if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
do concurrent (l = 1:3, m = 1:3)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
end do
do concurrent(l = 1:3, m = 1:3)
temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
@ -387,8 +387,8 @@ subroutine utilities_updateGamma(C)
call math_invert(A_inv, err, A)
temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* &
conjg(-xi1st(o,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset)
gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_complex(l,n)* &
conjg(-xi1st(o,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
end do
end if
end if
@ -405,7 +405,7 @@ end subroutine utilities_updateGamma
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorForward
tensorField_real(1:3,1:3,grid(1)+1:grid1Red*2,:,:) = 0.0_pReal
tensorField_real(1:3,1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
end subroutine utilities_FFTtensorForward
@ -429,7 +429,7 @@ end subroutine utilities_FFTtensorBackward
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarForward
scalarField_real(grid(1)+1:grid1Red*2,:,:) = 0.0_pReal
scalarField_real(cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
end subroutine utilities_FFTscalarForward
@ -454,7 +454,7 @@ end subroutine utilities_FFTscalarBackward
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorForward
vectorField_real(1:3,grid(1)+1:grid1Red*2,:,:) = 0.0_pReal
vectorField_real(1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
end subroutine utilities_FFTvectorForward
@ -493,8 +493,8 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium)
memoryEfficient: if (num%memory_efficient) then
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red
if (any([i,j,k+cells3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
do concurrent(l = 1:3, m = 1:3)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
end do
@ -519,7 +519,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
end if
end do; end do; end do
else memoryEfficient
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
do concurrent(l = 1:3, m = 1:3)
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k))
end do
@ -527,7 +527,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
end do; end do; end do
end if memoryEfficient
if (grid3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
end subroutine utilities_fourierGammaConvolution
@ -544,7 +544,7 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation
do k = 1, grid3; do j = 1, grid(2) ;do i = 1, grid1Red
do k = 1, cells3; do j = 1, cells(2) ;do i = 1, grid1Red
GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) &
/ (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal) &
* sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k))))
@ -571,7 +571,7 @@ real(pReal) function utilities_divergenceRMS()
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
utilities_divergenceRMS = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2)
do k = 1, cells3; do j = 1, cells(2)
do i = 2, grid1Red -1 ! Has somewhere a conj. complex counterpart. Therefore count it twice.
utilities_divergenceRMS = utilities_divergenceRMS &
+ 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k), & ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2, i.e. do not take square root and square again
@ -579,7 +579,7 @@ real(pReal) function utilities_divergenceRMS()
+sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),&
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2))
enddo
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1)
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if cells(1) /= 1)
+ sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
@ -589,7 +589,7 @@ real(pReal) function utilities_divergenceRMS()
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2)
enddo; enddo
if (grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
if (cells(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of cells(1) == 1
call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
@ -616,7 +616,7 @@ real(pReal) function utilities_curlRMS()
! calculating max curl criterion in Fourier space
utilities_curlRMS = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2);
do k = 1, cells3; do j = 1, cells(2);
do i = 2, grid1Red - 1
do l = 1, 3
curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) &
@ -638,7 +638,7 @@ real(pReal) function utilities_curlRMS()
-tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2))
enddo
utilities_curlRMS = utilities_curlRMS &
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1)
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (DC) does not have a conjugate complex counterpart (if cells(1) /= 1)
do l = 1, 3
curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3))
@ -648,13 +648,13 @@ real(pReal) function utilities_curlRMS()
-tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2))
enddo
utilities_curlRMS = utilities_curlRMS &
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1)
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (Nyquist) does not have a conjugate complex counterpart (if cells(1) /= 1)
enddo; enddo
call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
utilities_curlRMS = sqrt(utilities_curlRMS) * wgt
if (grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
if (cells(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of cells(1) == 1
end function utilities_curlRMS
@ -666,7 +666,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance
real(pReal), intent(in), dimension(3,3,3,3) :: C !< current average stiffness
type(rotation), intent(in) :: rot_BC !< rotation of load frame
type(tRotation), intent(in) :: rot_BC !< rotation of load frame
logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC
integer :: i, j
@ -736,7 +736,7 @@ subroutine utilities_fourierScalarGradient()
integer :: i, j, k
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) ! ToDo: no -conjg?
enddo; enddo; enddo
@ -750,7 +750,7 @@ subroutine utilities_fourierVectorDivergence()
integer :: i, j, k
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
scalarField_fourier(i,j,k) = sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k)))
enddo; enddo; enddo
@ -764,7 +764,7 @@ subroutine utilities_fourierVectorGradient()
integer :: i, j, k, m, n
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
do m = 1, 3; do n = 1, 3
tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k)
enddo; enddo
@ -780,7 +780,7 @@ subroutine utilities_fourierTensorDivergence()
integer :: i, j, k
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
vectorField_fourier(:,i,j,k) = matmul(tensorField_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k)))
enddo; enddo; enddo
@ -795,10 +795,10 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target
real(pReal), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress
real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: F !< deformation gradient target
real(pReal), intent(in) :: Delta_t !< loading time
type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame
type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame
integer :: i
@ -810,15 +810,15 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
print'(/,1x,a)', '... evaluating constitutive response ......................................'
flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field
call homogenization_mechanical_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field
call homogenization_mechanical_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) ! calculate P field
if (.not. terminallyIll) &
call homogenization_thermal_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3])
call homogenization_thermal_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
if (.not. terminallyIll) &
call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(grid(1:2))*grid3])
call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@ -833,7 +833,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
dPdF_norm_max = 0.0_pReal
dPdF_min = huge(1.0_pReal)
dPdF_norm_min = huge(1.0_pReal)
do i = 1, product(grid(1:2))*grid3
do i = 1, product(cells(1:2))*cells3
if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)
@ -878,16 +878,16 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
dt !< Delta_t between field0 and field
logical, intent(in) :: &
heterogeneous !< calculate field of rates
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
field0, & !< data of previous step
field !< data of current step
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: &
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_calculateRate
if (heterogeneous) then
utilities_calculateRate = (field-field0) / dt
else
utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3)
utilities_calculateRate = spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3)
endif
end function utilities_calculateRate
@ -901,12 +901,12 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim)
real(pReal), intent(in) :: &
Delta_t !< Delta_t of current step
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
field_lastInc, & !< initial field
rate !< rate by which to forward
real(pReal), intent(in), optional, dimension(3,3) :: &
aim !< average field value aim
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: &
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_forwardField
real(pReal), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
integer(MPI_INTEGER_KIND) :: err_MPI
@ -918,7 +918,7 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
fieldDiff = fieldDiff - aim
utilities_forwardField = utilities_forwardField - &
spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid3)
spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
endif
end function utilities_forwardField
@ -936,37 +936,37 @@ pure function utilities_getFreqDerivative(k_s)
select case (spectral_derivative_ID)
case (DERIVATIVE_CONTINUOUS_ID)
utilities_getFreqDerivative = cmplx(0.0_pReal, 2.0_pReal*PI*real(k_s,pReal)/geomSize,pReal)
utilities_getFreqDerivative = cmplx(0.0_pReal, TAU*real(k_s,pReal)/geomSize,pReal)
case (DERIVATIVE_CENTRAL_DIFF_ID)
utilities_getFreqDerivative = cmplx(0.0_pReal, sin(2.0_pReal*PI*real(k_s,pReal)/real(grid,pReal)), pReal)/ &
cmplx(2.0_pReal*geomSize/real(grid,pReal), 0.0_pReal, pReal)
utilities_getFreqDerivative = cmplx(0.0_pReal, sin(TAU*real(k_s,pReal)/real(cells,pReal)), pReal)/ &
cmplx(2.0_pReal*geomSize/real(cells,pReal), 0.0_pReal, pReal)
case (DERIVATIVE_FWBW_DIFF_ID)
utilities_getFreqDerivative(1) = &
cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) - 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(1)/real(grid(1),pReal), 0.0_pReal, pReal)
cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) - 1.0_pReal, &
sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(1)/real(cells(1),pReal), 0.0_pReal, pReal)
utilities_getFreqDerivative(2) = &
cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) - 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(2)/real(grid(2),pReal), 0.0_pReal, pReal)
cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) - 1.0_pReal, &
sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(2)/real(cells(2),pReal), 0.0_pReal, pReal)
utilities_getFreqDerivative(3) = &
cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* &
cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) - 1.0_pReal, &
sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(3)/real(grid(3),pReal), 0.0_pReal, pReal)
cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) + 1.0_pReal, &
sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* &
cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) - 1.0_pReal, &
sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ &
cmplx(4.0_pReal*geomSize(3)/real(cells(3),pReal), 0.0_pReal, pReal)
end select
end function utilities_getFreqDerivative
@ -979,10 +979,10 @@ end function utilities_getFreqDerivative
!--------------------------------------------------------------------------------------------------
subroutine utilities_updateCoords(F)
real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
real(pReal), dimension(3, grid(1),grid(2),grid3) :: IPcoords
real(pReal), dimension(3, grid(1),grid(2),grid3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI)
real(pReal), dimension(3, grid(1)+1,grid(2)+1,grid3+1) :: nodeCoords
real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: F
real(pReal), dimension(3, cells(1),cells(2),cells3) :: IPcoords
real(pReal), dimension(3, cells(1),cells(2),cells3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI)
real(pReal), dimension(3, cells(1)+1,cells(2)+1,cells3+1) :: nodeCoords
integer :: &
i,j,k,n, &
c
@ -1010,14 +1010,14 @@ subroutine utilities_updateCoords(F)
1, 1, 1, &
0, 1, 1 ], [3,8])
step = geomSize/real(grid, pReal)
step = geomSize/real(cells, pReal)
!--------------------------------------------------------------------------------------------------
! integration in Fourier space to get fluctuations of cell center discplacements
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F
call utilities_FFTtensorForward()
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any([i,j,k+grid3Offset] /= 1)) then
do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red
if (any([i,j,k+cells3Offset] /= 1)) then
vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) &
/ sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal)
else
@ -1029,25 +1029,25 @@ subroutine utilities_updateCoords(F)
!--------------------------------------------------------------------------------------------------
! average F
if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
if (cells3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
!--------------------------------------------------------------------------------------------------
! pad cell center fluctuations along z-direction (needed when running MPI simulation)
IPfluct_padded(1:3,1:grid(1),1:grid(2),2:grid3+1) = vectorField_real(1:3,1:grid(1),1:grid(2),1:grid3)
c = product(shape(IPfluct_padded(:,:,:,1))) !< amount of data to transfer
IPfluct_padded(1:3,1:cells(1),1:cells(2),2:cells3+1) = vectorField_real(1:3,1:cells(1),1:cells(2),1:cells3)
c = product(shape(IPfluct_padded(:,:,:,1))) !< amount of data to transfer
rank_t = modulo(worldrank+1_MPI_INTEGER_KIND,worldsize)
rank_b = modulo(worldrank-1_MPI_INTEGER_KIND,worldsize)
! send bottom layer to process below
call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(1),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI)
call MPI_Irecv(IPfluct_padded(:,:,:,cells3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
! send top layer to process above
call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI)
call MPI_Isend(IPfluct_padded(:,:,:,cells3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(4),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@ -1063,24 +1063,24 @@ subroutine utilities_updateCoords(F)
!--------------------------------------------------------------------------------------------------
! calculate nodal displacements
nodeCoords = 0.0_pReal
do k = 0,grid3; do j = 0,grid(2); do i = 0,grid(1)
nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+grid3Offset],pReal)))
do k = 0,cells3; do j = 0,cells(2); do i = 0,cells(1)
nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+cells3Offset],pReal)))
averageFluct: do n = 1,8
me = [i+neighbor(1,n),j+neighbor(2,n),k+neighbor(3,n)]
nodeCoords(1:3,i+1,j+1,k+1) = nodeCoords(1:3,i+1,j+1,k+1) &
+ IPfluct_padded(1:3,modulo(me(1)-1,grid(1))+1,modulo(me(2)-1,grid(2))+1,me(3)+1)*0.125_pReal
+ IPfluct_padded(1:3,modulo(me(1)-1,cells(1))+1,modulo(me(2)-1,cells(2))+1,me(3)+1)*0.125_pReal
enddo averageFluct
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! calculate cell center displacements
do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1)
do k = 1,cells3; do j = 1,cells(2); do i = 1,cells(1)
IPcoords(1:3,i,j,k) = vectorField_real(1:3,i,j,k) &
+ matmul(Favg,step*(real([i,j,k+grid3Offset],pReal)-0.5_pReal))
+ matmul(Favg,step*(real([i,j,k+cells3Offset],pReal)-0.5_pReal))
enddo; enddo; enddo
call discretization_setNodeCoords(reshape(NodeCoords,[3,(grid(1)+1)*(grid(2)+1)*(grid3+1)]))
call discretization_setIPcoords (reshape(IPcoords, [3,grid(1)*grid(2)*grid3]))
call discretization_setNodeCoords(reshape(NodeCoords,[3,(cells(1)+1)*(cells(2)+1)*(cells3+1)]))
call discretization_setIPcoords (reshape(IPcoords, [3,cells(1)*cells(2)*cells3]))
end subroutine utilities_updateCoords

4
src/homogenization.f90
View File

@ -414,7 +414,7 @@ subroutine homogenization_restartWrite(fileHandle)
groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_homogenization(ho))
call HDF5_write(homogState(ho)%state,groupHandle(2),'omega') ! ToDo: should be done by mech
call HDF5_write(homogState(ho)%state,groupHandle(2),'omega_mechanical') ! ToDo: should be done by mech
call HDF5_closeGroup(groupHandle(2))
@ -441,7 +441,7 @@ subroutine homogenization_restartRead(fileHandle)
groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_homogenization(ho))
call HDF5_read(homogState(ho)%state0,groupHandle(2),'omega') ! ToDo: should be done by mech
call HDF5_read(homogState(ho)%state0,groupHandle(2),'omega_mechanical') ! ToDo: should be done by mech
call HDF5_closeGroup(groupHandle(2))

18
src/lattice.f90
View File

@ -498,7 +498,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin
real(pReal), dimension(3,3,sum(Ntwin)):: coordinateSystem
type(rotation) :: R
type(tRotation) :: R
integer :: i
@ -538,7 +538,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66
real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S
type(rotation) :: R
type(tRotation) :: R
integer :: i
!--------------------------------------------------------------------------------------------------
@ -599,7 +599,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem !< coordinate system of slip system
real(pReal), dimension(3) :: direction, normal, np
type(rotation) :: R
type(tRotation) :: R
integer :: i
@ -1496,11 +1496,11 @@ function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_fcc,a_bcc) re
if (lattice_target == 'hP' .and. present(cOverA)) then
if (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal) &
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA=cOverA)
else if (lattice_target == 'cI' .and. present(a_fcc) .and. present(a_bcc)) then
if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) &
call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,a_fcc=a_fcc,a_bcc=a_bcc)
else
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
@ -1967,7 +1967,7 @@ end function buildCoordinateSystem
!--------------------------------------------------------------------------------------------------
subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
integer, dimension(:), intent(in) :: &
integer, dimension(:), intent(in) :: &
Ntrans
real(pReal), dimension(3,3,sum(Ntrans)), intent(out) :: &
Q, & !< Total rotation: Q = R*B
@ -1977,7 +1977,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
a_bcc, & !< lattice parameter a for bcc target lattice
a_fcc !< lattice parameter a for fcc parent lattice
type(rotation) :: &
type(tRotation) :: &
R, & !< Pitsch rotation
B !< Rotation of fcc to Bain coordinate system
real(pReal), dimension(3,3) :: &
@ -2079,10 +2079,10 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
Q(1:3,2,i) = y
Q(1:3,3,i) = z
S(1:3,1:3,i) = matmul(Q(1:3,1:3,i), matmul(matmul(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3 ! ToDo: This is of interest for the Schmid matrix only
enddo
end do
else
call IO_error(132,ext_msg='buildTransformationSystem')
endif
end if
end subroutine buildTransformationSystem

4
src/material.f90
View File

@ -18,7 +18,7 @@ module material
private
type :: tRotationContainer
type(Rotation), dimension(:), allocatable :: data
type(tRotation), dimension(:), allocatable :: data
end type
type :: tTensorContainer
real(pReal), dimension(:,:,:), allocatable :: data
@ -66,7 +66,7 @@ subroutine material_init(restart)
print'(/,1x,a)', '<<<+- material init -+>>>'; flush(IO_STDOUT)
call parse
call parse()
print'(/,1x,a)', 'parsed material.yaml'

55
src/math.f90
View File

@ -21,10 +21,11 @@ module math
config
#endif
real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter
real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian to degree
real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree to radian
complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi)
real(pReal), parameter :: &
PI = acos(-1.0_pReal), & !< ratio of a circle's circumference to its diameter
TAU = 2.0_pReal*PI, & !< ratio of a circle's circumference to its radius
INDEG = 360.0_pReal/TAU, & !< conversion from radian to degree
INRAD = TAU/360.0_pReal !< conversion from degree to radian
real(pReal), dimension(3,3), parameter :: &
math_I3 = reshape([&
@ -882,7 +883,7 @@ end function math_Voigt6to33_strain
!--------------------------------------------------------------------------------------------------
!> @brief Convert 3x3 tensor into 6 Voigt stress vector.
!> @brief Convert 3x3 stress tensor into 6 Voigt vector.
!--------------------------------------------------------------------------------------------------
pure function math_33toVoigt6_stress(sigma) result(sigma_tilde)
@ -897,7 +898,7 @@ end function math_33toVoigt6_stress
!--------------------------------------------------------------------------------------------------
!> @brief Convert 3x3 tensor into 6 Voigt strain vector.
!> @brief Convert 3x3 strain tensor into 6 Voigt vector.
!--------------------------------------------------------------------------------------------------
pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde)
@ -913,48 +914,48 @@ end function math_33toVoigt6_strain
!--------------------------------------------------------------------------------------------------
!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix.
!> @brief Convert 6x6 Voigt stiffness matrix into symmetric 3x3x3x3 tensor.
!--------------------------------------------------------------------------------------------------
pure function math_Voigt66to3333(m66)
pure function math_Voigt66to3333_stiffness(C_tilde) result(C)
real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333
real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix
real(pReal), dimension(3,3,3,3) :: C
real(pReal), dimension(6,6), intent(in) :: C_tilde
integer :: i,j
do i=1,6; do j=1,6
math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j)
math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j)
C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = C_tilde(i,j)
C(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = C_tilde(i,j)
C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = C_tilde(i,j)
C(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = C_tilde(i,j)
end do; end do
end function math_Voigt66to3333
end function math_Voigt66to3333_stiffness
!--------------------------------------------------------------------------------------------------
!> @brief Convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix.
!> @brief Convert 3x3x3x3 stiffness tensor into 6x6 Voigt matrix.
!--------------------------------------------------------------------------------------------------
pure function math_3333toVoigt66(m3333)
pure function math_3333toVoigt66_stiffness(C) result(C_tilde)
real(pReal), dimension(6,6) :: math_3333toVoigt66
real(pReal), dimension(3,3,3,3), intent(in) :: m3333 !< symmetric 3x3x3x3 matrix (no internal check)
real(pReal), dimension(6,6) :: C_tilde
real(pReal), dimension(3,3,3,3), intent(in) :: C
integer :: i,j
#ifndef __INTEL_COMPILER
do concurrent(i=1:6, j=1:6)
math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
end do
#else
do i=1,6; do j=1,6
math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
end do; end do
#endif
end function math_3333toVoigt66
end function math_3333toVoigt66_stiffness
!--------------------------------------------------------------------------------------------------
@ -984,7 +985,7 @@ impure elemental subroutine math_normal(x,mu,sigma)
end if
call random_number(rnd)
x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(2.0_pReal*PI*(1.0_pReal - rnd(2)))
x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(TAU*(1.0_pReal - rnd(2)))
end subroutine math_normal
@ -1088,7 +1089,7 @@ pure function math_rotationalPart(F) result(R)
if (dNeq0(x)) then
Phi = acos(math_clip((I_C(1)**3 -4.5_pReal*I_C(1)*I_C(2) +13.5_pReal*I_C(3))/x,-1.0_pReal,1.0_pReal))
lambda = I_C(1) +(2.0_pReal * sqrt(math_clip(I_C(1)**2-3.0_pReal*I_C(2),0.0_pReal))) &
*cos((Phi-2.0_pReal * PI*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal)
*cos((Phi-TAU*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal)
lambda = sqrt(math_clip(lambda,0.0_pReal)/3.0_pReal)
else
lambda = sqrt(I_C(1)/3.0_pReal)
@ -1154,8 +1155,8 @@ pure function math_eigvalsh33(m)
phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
math_eigvalsh33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
[cos( phi /3.0_pReal), &
cos((phi+2.0_pReal*PI)/3.0_pReal), &
cos((phi+4.0_pReal*PI)/3.0_pReal) &
cos((phi+TAU)/3.0_pReal), &
cos((phi+2.0_pReal*TAU)/3.0_pReal) &
] &
+ I(1)/3.0_pReal
endif
@ -1343,7 +1344,7 @@ subroutine selfTest
if (any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) &
error stop 'math_sym3333to66/math_66toSym3333'
if (any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) &
if (any(dNeq(math_3333toVoigt66_stiffness(math_Voigt66to3333_stiffness(t66)),t66,1.0e-15_pReal))) &
error stop 'math_3333toVoigt66/math_Voigt66to3333'
call random_number(v6)

4
src/mesh/DAMASK_mesh.f90
View File

@ -307,9 +307,11 @@ program DAMASK_mesh
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
end if
if (.not. cutBack .and. worldrank == 0) &
if (.not. cutBack .and. worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
flush(statUnit)
endif
end do subStepLooping
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc

20
src/parallelization.f90
View File

@ -52,13 +52,13 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine parallelization_init
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize, version, subversion, devNull
character(len=4) :: rank_str
character(len=MPI_MAX_LIBRARY_VERSION_STRING) :: MPI_library_version
!$ integer :: got_env, threadLevel
!$ integer(pI32) :: OMP_NUM_THREADS
!$ character(len=6) NumThreadsString
PetscErrorCode :: err_PETSc
#ifdef _OPENMP
! If openMP is enabled, check if the MPI libary supports it and initialize accordingly.
@ -86,12 +86,22 @@ subroutine parallelization_init
if (err_MPI /= 0_MPI_INTEGER_KIND) &
error stop 'Could not determine worldrank'
if (worldrank == 0) print'(/,1x,a)', '<<<+- parallelization init -+>>>'
if (worldrank == 0) then
print'(/,1x,a)', '<<<+- parallelization init -+>>>'
call MPI_Get_library_version(MPI_library_version,devNull,err_MPI)
print'(/,1x,a)', trim(MPI_library_version)
call MPI_Get_version(version,subversion,err_MPI)
print'(1x,a,i0,a,i0)', 'MPI standard: ',version,'.',subversion
#ifdef _OPENMP
print'(1x,a,i0)', 'OpenMP version: ',openmp_version
#endif
end if
call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) &
error stop 'Could not determine worldsize'
if (worldrank == 0) print'(/,1x,a,i3)', 'MPI processes: ',worldsize
if (worldrank == 0) print'(/,1x,a,i0)', 'MPI processes: ',worldsize
call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) &
@ -128,7 +138,7 @@ subroutine parallelization_init
!$ OMP_NUM_THREADS = 4_pI32
!$ endif
!$ endif
!$ print'(1x,a,1x,i2)', 'OMP_NUM_THREADS:',OMP_NUM_THREADS
!$ print'(1x,a,i0)', 'OMP_NUM_THREADS: ',OMP_NUM_THREADS
!$ call omp_set_num_threads(OMP_NUM_THREADS)
end subroutine parallelization_init

50
src/phase.f90
View File

@ -8,6 +8,7 @@ module phase
use constants
use math
use rotations
use polynomials
use IO
use config
use material
@ -123,11 +124,20 @@ module phase
integer, intent(in) :: ph
end subroutine mechanical_restartWrite
module subroutine thermal_restartWrite(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph
end subroutine thermal_restartWrite
module subroutine mechanical_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph
end subroutine mechanical_restartRead
module subroutine thermal_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph
end subroutine thermal_restartRead
module function mechanical_S(ph,en) result(S)
integer, intent(in) :: ph,en
@ -313,7 +323,6 @@ module phase
phase_restore, &
plastic_nonlocal_updateCompatibility, &
converged, &
crystallite_init, &
phase_mechanical_constitutive, &
phase_thermal_constitutive, &
phase_damage_constitutive, &
@ -391,6 +400,8 @@ subroutine phase_init
call damage_init
call thermal_init(phases)
call crystallite_init()
end subroutine phase_init
@ -398,7 +409,7 @@ end subroutine phase_init
!> @brief Allocate the components of the state structure for a given phase
!--------------------------------------------------------------------------------------------------
subroutine phase_allocateState(state, &
NEntries,sizeState,sizeDotState,sizeDeltaState)
NEntries,sizeState,sizeDotState,sizeDeltaState,offsetDeltaState)
class(tState), intent(inout) :: &
state
@ -407,12 +418,17 @@ subroutine phase_allocateState(state, &
sizeState, &
sizeDotState, &
sizeDeltaState
integer, intent(in), optional :: &
offsetDeltaState
state%sizeState = sizeState
state%sizeDotState = sizeDotState
state%sizeDeltaState = sizeDeltaState
state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
if (present(offsetDeltaState)) then
state%offsetDeltaState = offsetDeltaState ! ToDo: this is a fix for broken nonlocal
else
state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
end if
allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,NEntries), source=0.0_pReal)
@ -421,7 +437,8 @@ subroutine phase_allocateState(state, &
allocate(state%dotState (sizeDotState,NEntries), source=0.0_pReal)
allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal)
state%deltaState2 => state%state(state%offsetDeltaState+1: &
state%offsetDeltaState+state%sizeDeltaState,:)
end subroutine phase_allocateState
@ -486,22 +503,13 @@ subroutine crystallite_init()
ce, &
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el, & !< counter in element loop
cMax, & !< maximum number of integration point components
iMax, & !< maximum number of integration points
eMax !< maximum number of elements
el !< counter in element loop
class(tNode), pointer :: &
num_crystallite, &
phases
print'(/,1x,a)', '<<<+- crystallite init -+>>>'
cMax = homogenization_maxNconstituents
iMax = discretization_nIPs
eMax = discretization_Nelems
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal)
@ -535,15 +543,9 @@ subroutine crystallite_init()
phases => config_material%get('phase')
print'(/,a42,1x,i10)', ' # of elements: ', eMax
print'( a42,1x,i10)', ' # of integration points/element: ', iMax
print'( a42,1x,i10)', 'max # of constituents/integration point: ', cMax
flush(IO_STDOUT)
!$OMP PARALLEL DO PRIVATE(ce)
do el = 1, eMax
do ip = 1, iMax
do el = 1, discretization_Nelems
do ip = 1, discretization_nIPs
ce = (el-1)*discretization_nIPs + ip
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
call crystallite_orientations(co,ip,el)
@ -640,6 +642,7 @@ subroutine phase_restartWrite(fileHandle)
groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_phase(ph))
call mechanical_restartWrite(groupHandle(2),ph)
call thermal_restartWrite(groupHandle(2),ph)
call HDF5_closeGroup(groupHandle(2))
@ -668,6 +671,7 @@ subroutine phase_restartRead(fileHandle)
groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_phase(ph))
call mechanical_restartRead(groupHandle(2),ph)
call thermal_restartRead(groupHandle(2),ph)
call HDF5_closeGroup(groupHandle(2))

152
src/phase_mechanical.f90
View File

@ -40,8 +40,6 @@ submodule(phase) mechanical
integer(kind(PLASTIC_undefined_ID)), dimension(:), allocatable :: &
phase_plasticity !< plasticity of each phase
integer :: phase_plasticity_maxSizeDotState
interface
module subroutine eigen_init(phases)
@ -81,16 +79,17 @@ submodule(phase) mechanical
en
end subroutine plastic_isotropic_LiAndItsTangent
module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
co, & !< constituent
ip, & !< integration point
el, & !< element
ph, &
en
real(pReal), intent(in) :: &
subdt !< timestep
logical :: broken
subdt !< timestep
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function plastic_dotState
module function plastic_deltaState(ph, en) result(broken)
@ -296,8 +295,6 @@ module subroutine mechanical_init(phases)
do ph = 1,phases%length
plasticState(ph)%state0 = plasticState(ph)%state
enddo
phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -588,9 +585,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: &
broken
@ -601,43 +598,43 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
sizeDotState
real(pReal) :: &
zeta
real(pReal), dimension(phase_plasticity_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(phase_plasticity_maxSizeDotState,2) :: &
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
r, & ! state residuum
dotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: &
dotState_last
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState * Delta_t
iteration: do NiterationState = 1, num%nState
dotState(1:sizeDotState,2) = merge(dotState(1:sizeDotState,1),0.0, nIterationState > 1)
dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,en)
dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1)
dotState_last(1:sizeDotState,1) = dotState
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
if(broken) exit iteration
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) exit iteration
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) exit iteration
zeta = damper(plasticState(ph)%dotState(:,en),dotState(1:sizeDotState,1),&
dotState(1:sizeDotState,2))
plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) * zeta &
+ dotState(1:sizeDotState,1) * (1.0_pReal - zeta)
r(1:sizeDotState) = plasticState(ph)%state(1:sizeDotState,en) &
- subState0 &
- plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) &
- r(1:sizeDotState)
if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then
zeta = damper(dotState,dotState_last(1:sizeDotState,1),dotState_last(1:sizeDotState,2))
dotState = dotState * zeta &
+ dotState_last(1:sizeDotState,1) * (1.0_pReal - zeta)
r = plasticState(ph)%state(1:sizeDotState,en) &
- subState0 &
- dotState * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) - r
if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then
broken = plastic_deltaState(ph,en)
exit iteration
endif
@ -652,19 +649,20 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
!--------------------------------------------------------------------------------------------------
real(pReal) pure function damper(omega_0,omega_1,omega_2)
real(pReal), dimension(:), intent(in) :: &
omega_0, omega_1, omega_2
real(pReal), dimension(:), intent(in) :: &
omega_0, omega_1, omega_2
real(pReal) :: dot_prod12, dot_prod22
real(pReal) :: dot_prod12, dot_prod22
dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2)
dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2)
if (min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) then
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
else
damper = 1.0_pReal
endif
dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2)
dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2)
if (min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) then
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
else
damper = 1.0_pReal
endif
end function damper
@ -686,6 +684,8 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
logical :: &
broken
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
dotState
integer :: &
ph, &
en, &
@ -694,13 +694,14 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t
+ dotState * Delta_t
broken = plastic_deltaState(ph,en)
if(broken) return
@ -729,20 +730,23 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
ph, &
en, &
sizeDotState
real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
r, &
dotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,en) * 0.5_pReal * Delta_t
r = - dotState * 0.5_pReal * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotstate(1:sizeDotState,en) * Delta_t
+ dotState * Delta_t
broken = plastic_deltaState(ph,en)
if(broken) return
@ -750,10 +754,10 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
if(broken) return
broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,en) * Delta_t, &
broken = .not. converged(r + 0.5_pReal * dotState * Delta_t, &
plasticState(ph)%state(1:sizeDotState,en), &
plasticState(ph)%atol(1:sizeDotState))
@ -847,44 +851,48 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
ph, &
en, &
sizeDotState
real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
dotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: &
plastic_RKdotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
broken = plastic_dotState(Delta_t,co,ip,el,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
do stage = 1, size(A,1)
plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,en)
plasticState(ph)%dotState(:,en) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
plastic_RKdotState(1:sizeDotState,stage) = dotState
dotState = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
do n = 2, stage
plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) &
+ A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
dotState = dotState &
+ A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
enddo
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t
+ dotState * Delta_t
broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el)
broken = integrateStress(F_0 + (F-F_0) * Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage),co,ip,el)
if(broken) exit
broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,en)
if(broken) exit
dotState = plastic_dotState(Delta_t*C(stage), co,ip,el,ph,en)
if (any(IEEE_is_NaN(dotState))) exit
enddo
if(broken) return
plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,en)
plasticState(ph)%dotState(:,en) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B)
plastic_RKdotState(1:sizeDotState,size(B)) = dotState
dotState = matmul(plastic_RKdotState,B)
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t
+ dotState * Delta_t
if(present(DB)) &
broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, &
@ -958,7 +966,7 @@ subroutine crystallite_results(group,ph)
!--------------------------------------------------------------------------------------------------
function to_quaternion(dataset)
type(rotation), dimension(:), intent(in) :: dataset
type(tRotation), dimension(:), intent(in) :: dataset
real(pReal), dimension(4,size(dataset,1)) :: to_quaternion
integer :: i
@ -1248,7 +1256,7 @@ module subroutine mechanical_restartWrite(groupHandle,ph)
integer, intent(in) :: ph
call HDF5_write(plasticState(ph)%state,groupHandle,'omega')
call HDF5_write(plasticState(ph)%state,groupHandle,'omega_plastic')
call HDF5_write(phase_mechanical_Fi(ph)%data,groupHandle,'F_i')
call HDF5_write(phase_mechanical_Li(ph)%data,groupHandle,'L_i')
call HDF5_write(phase_mechanical_Lp(ph)%data,groupHandle,'L_p')
@ -1265,7 +1273,7 @@ module subroutine mechanical_restartRead(groupHandle,ph)
integer, intent(in) :: ph
call HDF5_read(plasticState(ph)%state0,groupHandle,'omega')
call HDF5_read(plasticState(ph)%state0,groupHandle,'omega_plastic')
call HDF5_read(phase_mechanical_Fi0(ph)%data,groupHandle,'F_i')
call HDF5_read(phase_mechanical_Li0(ph)%data,groupHandle,'L_i')
call HDF5_read(phase_mechanical_Lp0(ph)%data,groupHandle,'L_p')

47
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -8,10 +8,9 @@ submodule(phase:eigen) thermalexpansion
integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance
type :: tParameters
real(pReal) :: &
T_ref
real(pReal), dimension(3,3,3) :: &
A = 0.0_pReal
type(tPolynomial) :: &
A_11, &
A_33
end type tParameters
type(tParameters), dimension(:), allocatable :: param
@ -28,13 +27,13 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstances,p,i,k
integer :: Ninstances, p, k
class(tNode), pointer :: &
phases, &
phase, &
mech, &
kinematics, &
kinematic_type
myConfig
print'(/,1x,a)', '<<<+- phase:mechanical:eigen:thermalexpansion init -+>>>'
@ -56,21 +55,13 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
do k = 1, kinematics%length
if (myKinematics(k,p)) then
associate(prm => param(kinematics_thermal_expansion_instance(p)))
kinematic_type => kinematics%get(k)
prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=T_ROOM)
myConfig => kinematics%get(k)
prm%A_11 = polynomial(myConfig%asDict(),'A_11','T')
if (any(phase_lattice(p) == ['hP','tI'])) &
prm%A_33 = polynomial(myConfig%asDict(),'A_33','T')
prm%A(1,1,1) = kinematic_type%get_asFloat('A_11')
prm%A(1,1,2) = kinematic_type%get_asFloat('A_11,T', defaultVal=0.0_pReal)
prm%A(1,1,3) = kinematic_type%get_asFloat('A_11,T^2',defaultVal=0.0_pReal)
if (any(phase_lattice(p) == ['hP','tI'])) then
prm%A(3,3,1) = kinematic_type%get_asFloat('A_33')
prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T', defaultVal=0.0_pReal)
prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal)
end if
do i=1, size(prm%A,3)
prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p))
end do
end associate
end if
end do
@ -91,22 +82,20 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
real(pReal) :: T, dot_T
real(pReal), dimension(3,3) :: A
T = thermal_T(ph,me)
dot_T = thermal_dot_T(ph,me)
associate(prm => param(kinematics_thermal_expansion_instance(ph)))
Li = dot_T * ( &
prm%A(1:3,1:3,1) & ! constant coefficient
+ prm%A(1:3,1:3,2)*(T - prm%T_ref) & ! linear coefficient
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
) / &
(1.0_pReal &
+ prm%A(1:3,1:3,1)*(T - prm%T_ref) / 1.0_pReal &
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2.0_pReal &
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3.0_pReal &
)
A = 0.0_pReal
A(1,1) = prm%A_11%at(T)
if (any(phase_lattice(ph) == ['hP','tI'])) A(3,3) = prm%A_33%at(T)
A = lattice_symmetrize_33(A,phase_lattice(ph))
Li = dot_T * A
end associate
dLi_dTstar = 0.0_pReal

82
src/phase_mechanical_elastic.f90
View File

@ -1,15 +1,13 @@
submodule(phase:mechanical) elastic
type :: tParameters
real(pReal),dimension(3) :: &
C_11 = 0.0_pReal, &
C_12 = 0.0_pReal, &
C_13 = 0.0_pReal, &
C_33 = 0.0_pReal, &
C_44 = 0.0_pReal, &
C_66 = 0.0_pReal
real(pReal) :: &
T_ref
type(tPolynomial) :: &
C_11, &
C_12, &
C_13, &
C_33, &
C_44, &
C_66
end type tParameters
type(tParameters), allocatable, dimension(:) :: param
@ -47,35 +45,17 @@ module subroutine elastic_init(phases)
associate(prm => param(ph))
prm%T_ref = elastic%get_asFloat('T_ref', defaultVal=T_ROOM)
prm%C_11(1) = elastic%get_asFloat('C_11')
prm%C_11(2) = elastic%get_asFloat('C_11,T', defaultVal=0.0_pReal)
prm%C_11(3) = elastic%get_asFloat('C_11,T^2',defaultVal=0.0_pReal)
prm%C_12(1) = elastic%get_asFloat('C_12')
prm%C_12(2) = elastic%get_asFloat('C_12,T', defaultVal=0.0_pReal)
prm%C_12(3) = elastic%get_asFloat('C_12,T^2',defaultVal=0.0_pReal)
prm%C_44(1) = elastic%get_asFloat('C_44')
prm%C_44(2) = elastic%get_asFloat('C_44,T', defaultVal=0.0_pReal)
prm%C_44(3) = elastic%get_asFloat('C_44,T^2',defaultVal=0.0_pReal)
prm%C_11 = polynomial(elastic%asDict(),'C_11','T')
prm%C_12 = polynomial(elastic%asDict(),'C_12','T')
prm%C_44 = polynomial(elastic%asDict(),'C_44','T')
if (any(phase_lattice(ph) == ['hP','tI'])) then
prm%C_13(1) = elastic%get_asFloat('C_13')
prm%C_13(2) = elastic%get_asFloat('C_13,T', defaultVal=0.0_pReal)
prm%C_13(3) = elastic%get_asFloat('C_13,T^2',defaultVal=0.0_pReal)
prm%C_33(1) = elastic%get_asFloat('C_33')
prm%C_33(2) = elastic%get_asFloat('C_33,T', defaultVal=0.0_pReal)
prm%C_33(3) = elastic%get_asFloat('C_33,T^2',defaultVal=0.0_pReal)
prm%C_13 = polynomial(elastic%asDict(),'C_13','T')
prm%C_33 = polynomial(elastic%asDict(),'C_33','T')
end if
if (phase_lattice(ph) == 'tI') then
prm%C_66(1) = elastic%get_asFloat('C_66')
prm%C_66(2) = elastic%get_asFloat('C_66,T', defaultVal=0.0_pReal)
prm%C_66(3) = elastic%get_asFloat('C_66,T^2',defaultVal=0.0_pReal)
end if
if (phase_lattice(ph) == 'tI') &
prm%C_66 = polynomial(elastic%asDict(),'C_66','T')
end associate
end do
@ -97,38 +77,20 @@ pure module function elastic_C66(ph,en) result(C66)
associate(prm => param(ph))
C66 = 0.0_pReal
T = thermal_T(ph,en)
C66(1,1) = prm%C_11(1) &
+ prm%C_11(2)*(T - prm%T_ref) &
+ prm%C_11(3)*(T - prm%T_ref)**2
C66(1,2) = prm%C_12(1) &
+ prm%C_12(2)*(T - prm%T_ref) &
+ prm%C_12(3)*(T - prm%T_ref)**2
C66(4,4) = prm%C_44(1) &
+ prm%C_44(2)*(T - prm%T_ref) &
+ prm%C_44(3)*(T - prm%T_ref)**2
C66(1,1) = prm%C_11%at(T)
C66(1,2) = prm%C_12%at(T)
C66(4,4) = prm%C_44%at(T)
if (any(phase_lattice(ph) == ['hP','tI'])) then
C66(1,3) = prm%C_13(1) &
+ prm%C_13(2)*(T - prm%T_ref) &
+ prm%C_13(3)*(T - prm%T_ref)**2
C66(3,3) = prm%C_33(1) &
+ prm%C_33(2)*(T - prm%T_ref) &
+ prm%C_33(3)*(T - prm%T_ref)**2
C66(1,3) = prm%C_13%at(T)
C66(3,3) = prm%C_33%at(T)
end if
if (phase_lattice(ph) == 'tI') then
C66(6,6) = prm%C_66(1) &
+ prm%C_66(2)*(T - prm%T_ref) &
+ prm%C_66(3)*(T - prm%T_ref)**2
end if
if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66%at(T)
C66 = lattice_symmetrize_C66(C66,phase_lattice(ph))
@ -200,7 +162,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
C66 = phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)
C = math_Voigt66to3333(C66)
C = math_Voigt66to3333_stiffness(C66)
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_Voigt6to33_stress(matmul(C66,math_33toVoigt6_strain(matmul(matmul(transpose(Fi),E),Fi))))!< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration

59
src/phase_mechanical_plastic.f90
View File

@ -110,29 +110,35 @@ submodule(phase:mechanical) plastic
end subroutine nonlocal_LpAndItsTangent
module subroutine isotropic_dotState(Mp,ph,en)
module function isotropic_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine isotropic_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function isotropic_dotState
module subroutine phenopowerlaw_dotState(Mp,ph,en)
module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine phenopowerlaw_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function phenopowerlaw_dotState
module subroutine plastic_kinehardening_dotState(Mp,ph,en)
module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine plastic_kinehardening_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function plastic_kinehardening_dotState
module subroutine dislotwin_dotState(Mp,T,ph,en)
real(pReal), dimension(3,3), intent(in) :: &
@ -144,21 +150,20 @@ submodule(phase:mechanical) plastic
en
end subroutine dislotwin_dotState
module subroutine dislotungsten_dotState(Mp,T,ph,en)
module function dislotungsten_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
ph, &
en
end subroutine dislotungsten_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function dislotungsten_dotState
module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,en,ip,el)
module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
real(pReal), intent(in) :: &
Temperature, & !< temperature
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
@ -283,7 +288,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
do i=1,3; do j=1,3
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
end do; end do
end if
@ -294,7 +299,7 @@ end subroutine plastic_LpAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
integer, intent(in) :: &
co, & !< component-ID of integration point
@ -306,7 +311,8 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
subdt !< timestep
real(pReal), dimension(3,3) :: &
Mp
logical :: broken
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then
@ -316,27 +322,28 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
plasticType: select case (phase_plasticity(ph))
case (PLASTIC_ISOTROPIC_ID) plasticType
call isotropic_dotState(Mp,ph,en)
dotState = isotropic_dotState(Mp,ph,en)
case (PLASTIC_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_dotState(Mp,ph,en)
dotState = phenopowerlaw_dotState(Mp,ph,en)
case (PLASTIC_KINEHARDENING_ID) plasticType
call plastic_kinehardening_dotState(Mp,ph,en)
dotState = plastic_kinehardening_dotState(Mp,ph,en)
case (PLASTIC_DISLOTWIN_ID) plasticType
call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en)
dotState = plasticState(ph)%dotState(:,en)
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en)
dotState = dislotungsten_dotState(Mp,ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType
call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el)
call nonlocal_dotState(Mp,subdt,ph,en,ip,el)
dotState = plasticState(ph)%dotState(:,en)
end select plasticType
end if
broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,en)))
end function plastic_dotState
@ -391,6 +398,7 @@ module function plastic_deltaState(ph, en) result(broken)
myOffset, &
mySize
broken = .false.
select case (phase_plasticity(ph))
@ -412,10 +420,9 @@ module function plastic_deltaState(ph, en) result(broken)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
if (.not. broken) then
myOffset = plasticState(ph)%offsetDeltaState
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = &
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en)
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) &
+ plasticState(ph)%deltaState(1:mySize,en)
end if
end select

59
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -43,6 +43,13 @@ submodule(phase:plastic) dislotungsten
systems_sl
end type tParameters !< container type for internal constitutive parameters
type :: tIndexDotState
integer, dimension(2) :: &
rho_mob, &
rho_dip, &
gamma_sl
end type tIndexDotState
type :: tDislotungstenState
real(pReal), dimension(:,:), pointer :: &
rho_mob, &
@ -58,10 +65,9 @@ submodule(phase:plastic) dislotungsten
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tDisloTungstenState), allocatable, dimension(:) :: &
dotState, &
state
type(tParameters), allocatable, dimension(:) :: param
type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tDisloTungstenState), allocatable, dimension(:) :: state
type(tDisloTungstenDependentState), allocatable, dimension(:) :: dependentState
contains
@ -103,18 +109,17 @@ module function plastic_dislotungsten_init() result(myPlasticity)
print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002'
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(dependentState(phases%length))
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
@ -214,28 +219,29 @@ module function plastic_dislotungsten_init() result(myPlasticity)
sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
idx_dot%rho_mob = [startIndex,endIndex]
stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob = spread(rho_mob_0,2,Nmembers)
dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
idx_dot%rho_dip = [startIndex,endIndex]
stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip = spread(rho_dip_0,2,Nmembers)
dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
idx_dot%gamma_sl = [startIndex,endIndex]
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
@ -300,15 +306,15 @@ end subroutine dislotungsten_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine dislotungsten_dotState(Mp,T,ph,en)
module function dislotungsten_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T !< temperature
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_pos, dot_gamma_neg,&
@ -319,17 +325,22 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot_rho_dip_climb, &
d_hat
real(pReal) :: &
mu
mu, T
associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
dot_rho_mob => dotState(indexDotState(ph)%rho_mob(1):indexDotState(ph)%rho_mob(2)), &
dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), &
dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)))
mu = elastic_mu(ph,en)
T = thermal_T(ph,en)
call kinetics(Mp,T,ph,en,&
dot_gamma_pos,dot_gamma_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg)
dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
dot_gamma_sl = abs(dot_gamma_pos+dot_gamma_neg)
where(dEq0((tau_pos+tau_neg)*0.5_pReal))
dot_rho_dip_formation = 0.0_pReal
@ -338,24 +349,24 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
d_hat = math_clip(3.0_pReal*mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), &
prm%d_caron, & ! lower limit
dst%Lambda_sl(:,en)) ! upper limit
dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, &
dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot_gamma_sl/prm%b_sl, &
0.0_pReal, &
prm%dipoleformation)
v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(2.0_pReal*PI*K_B*T)) &
v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(TAU*K_B*T)) &
* (1.0_pReal/(d_hat+prm%d_caron))
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency?
end where
dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication
dot_rho_mob = dot_gamma_sl/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication
- dot_rho_dip_formation &
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges
dot%rho_dip(:,en) = dot_rho_dip_formation &
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot_gamma_sl ! Spontaneous annihilation of 2 edges
dot_rho_dip = dot_rho_dip_formation &
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot_gamma_sl & ! Spontaneous annihilation of an edge with a dipole
- dot_rho_dip_climb
end associate
end subroutine dislotungsten_dotState
end function dislotungsten_dotState
!--------------------------------------------------------------------------------------------------

44
src/phase_mechanical_plastic_isotropic.f90
View File

@ -37,9 +37,7 @@ submodule(phase:plastic) isotropic
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tIsotropicState), allocatable, dimension(:) :: &
dotState, &
state
type(tIsotropicState), allocatable, dimension(:) :: state
contains
@ -77,16 +75,15 @@ module function plastic_isotropic_init() result(myPlasticity)
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
associate(prm => param(ph), stt => state(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
#if defined (__GFORTRAN__)
prm%output = output_as1dString(pl)
@ -125,12 +122,12 @@ module function plastic_isotropic_init() result(myPlasticity)
sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
stt%xi => plasticState(ph)%state (1,:)
stt%xi = xi_0
dot%xi => plasticState(ph)%dotState(1,:)
stt%xi => plasticState(ph)%state(1,:)
stt%xi = xi_0
plasticState(ph)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if (plasticState(ph)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
@ -178,7 +175,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
if (norm_Mp_dev > 0.0_pReal) then
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en)))**prm%n
Lp = dot_gamma * Mp_dev/norm_Mp_dev
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -242,27 +239,26 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine isotropic_dotState(Mp,ph,en)
module function isotropic_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal) :: &
dot_gamma, & !< strainrate
xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(ph), stt => state(ph), &
dot => dotState(ph))
associate(prm => param(ph), stt => state(ph), dot_xi => dotState(1))
if (prm%dilatation) then
norm_Mp = sqrt(math_tensordot(Mp,Mp))
else
norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp)))
end if
norm_Mp = merge(sqrt(math_tensordot(Mp,Mp)), &
sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))), &
prm%dilatation)
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n
@ -274,16 +270,16 @@ module subroutine isotropic_dotState(Mp,ph,en)
+ asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) &
/ prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
end if
dot%xi(en) = dot_gamma &
* ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star)
dot_xi = dot_gamma &
* ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star)
else
dot%xi(en) = 0.0_pReal
dot_xi = 0.0_pReal
end if
end associate
end subroutine isotropic_dotState
end function isotropic_dotState
!--------------------------------------------------------------------------------------------------

56
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -34,6 +34,13 @@ submodule(phase:plastic) kinehardening
systems_sl
end type tParameters
type :: tIndexDotState
integer, dimension(2) :: &
xi, &
chi, &
gamma
end type tIndexDotState
type :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: &
xi, & !< resistance against plastic slip
@ -47,10 +54,8 @@ submodule(phase:plastic) kinehardening
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tKinehardeningState), allocatable, dimension(:) :: &
dotState, &
deltaState, &
state
type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tKinehardeningState), allocatable, dimension(:) :: state, deltaState
contains
@ -91,19 +96,20 @@ module function plastic_kinehardening_init() result(myPlasticity)
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(deltaState(phases%length))
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph))
associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
#if defined (__GFORTRAN__)
prm%output = output_as1dString(pl)
@ -173,27 +179,28 @@ module function plastic_kinehardening_init() result(myPlasticity)
sizeState = sizeDotState + sizeDeltaState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi => plasticState(ph)%state (startIndex:endIndex,:)
idx_dot%xi = [startIndex,endIndex]
stt%xi => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi = spread(xi_0, 2, Nmembers)
dot%xi => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%chi => plasticState(ph)%state (startIndex:endIndex,:)
dot%chi => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%chi = [startIndex,endIndex]
stt%chi => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%gamma = [startIndex,endIndex]
stt%gamma => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
@ -270,13 +277,15 @@ end subroutine kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_dotState(Mp,ph,en)
module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal) :: &
sumGamma
@ -284,29 +293,32 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en)
dot_gamma_pos,dot_gamma_neg
associate(prm => param(ph), stt => state(ph),dot => dotState(ph))
associate(prm => param(ph), stt => state(ph), &
dot_xi => dotState(IndexDotState(ph)%xi(1):IndexDotState(ph)%xi(2)),&
dot_chi => dotState(IndexDotState(ph)%chi(1):IndexDotState(ph)%chi(2)),&
dot_gamma => dotState(IndexDotState(ph)%gamma(1):IndexDotState(ph)%gamma(2)))
call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
dot%gamma(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
dot_gamma = abs(dot_gamma_pos+dot_gamma_neg)
sumGamma = sum(stt%gamma(:,en))
dot%xi(:,en) = matmul(prm%h_sl_sl,dot%gamma(:,en)) &
dot_xi = matmul(prm%h_sl_sl,dot_gamma) &
* ( prm%h_inf_f &
+ (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
* exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
)
dot%chi(:,en) = stt%sgn_gamma(:,en)*dot%gamma(:,en) * &
( prm%h_inf_b + &
(prm%h_0_b - prm%h_inf_b &
dot_chi = stt%sgn_gamma(:,en)*dot_gamma &
* ( prm%h_inf_b &
+ (prm%h_0_b - prm%h_inf_b &
+ prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi_0(:,en))*(stt%gamma(:,en)-stt%gamma_0(:,en))&
) *exp(-(stt%gamma(:,en)-stt%gamma_0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi_0(:,en))) &
)
end associate
end subroutine plastic_kinehardening_dotState
end function plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------

22
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -227,8 +227,8 @@ module function plastic_nonlocal_init() result(myPlasticity)
st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
plasticState(ph)%nonlocal = pl%get_asBool('flux',defaultVal=.True.)
#if defined (__GFORTRAN__)
@ -403,7 +403,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention
allocate(geom(ph)%V_0(Nmembers))
call storeGeometry(ph)
@ -411,9 +411,6 @@ module function plastic_nonlocal_init() result(myPlasticity)
if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
call IO_error(212,ext_msg='IPneighborhood does not exist')
plasticState(ph)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
st0%rho => plasticState(ph)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho => plasticState(ph)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho => plasticState(ph)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
@ -941,13 +938,12 @@ end subroutine plastic_nonlocal_deltaState
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
module subroutine nonlocal_dotState(Mp,timestep, &
ph,en,ip,el)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
real(pReal), intent(in) :: &
Temperature, & !< temperature
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
@ -984,7 +980,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
real(pReal) :: &
D_SD, &
mu, &
nu
nu, Temperature
if (timestep <= 0.0_pReal) then
plasticState(ph)%dotState = 0.0_pReal
@ -995,6 +991,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en)
Temperature = thermal_T(ph,en)
tau = 0.0_pReal
dot_gamma = 0.0_pReal
@ -1195,7 +1192,6 @@ function rhoDotFlux(timestep,ph,en,ip,el)
associate(prm => param(ph), &
dst => dependentState(ph), &
dot => dotState(ph), &
stt => state(ph))
ns = prm%sum_N_sl
@ -1206,7 +1202,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
rho0 = getRho0(ph,en)
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here
dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
@ -1217,7 +1213,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
if (plasticState(ph)%nonlocal) then
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
.and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG
@ -1395,7 +1391,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
nThresholdValues
logical, dimension(param(ph)%sum_N_sl) :: &
belowThreshold
type(rotation) :: mis
type(tRotation) :: mis
associate(prm => param(ph))

70
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -47,6 +47,14 @@ submodule(phase:plastic) phenopowerlaw
systems_tw
end type tParameters
type :: tIndexDotState
integer, dimension(2) :: &
xi_sl, &
xi_tw, &
gamma_sl, &
gamma_tw
end type tIndexDotState
type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
xi_sl, &
@ -56,11 +64,10 @@ submodule(phase:plastic) phenopowerlaw
end type tPhenopowerlawState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
! containers for parameters, dot state index, and state
type(tParameters), allocatable, dimension(:) :: param
type(tPhenopowerlawState), allocatable, dimension(:) :: &
dotState, &
state
type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tPhenopowerlawState), allocatable, dimension(:) :: state
contains
@ -101,17 +108,18 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
associate(prm => param(ph), stt => state(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
!--------------------------------------------------------------------------------------------------
! slip related parameters
@ -224,37 +232,37 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi_sl => plasticState(ph)%state (startIndex:endIndex,:)
idx_dot%xi_sl = [startIndex,endIndex]
stt%xi_sl => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_sl = spread(xi_0_sl, 2, Nmembers)
dot%xi_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%xi_tw => plasticState(ph)%state (startIndex:endIndex,:)
idx_dot%xi_tw = [startIndex,endIndex]
stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_tw = spread(xi_0_tw, 2, Nmembers)
dot%xi_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%gamma_sl = [startIndex,endIndex]
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%gamma_tw => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%gamma_tw = [startIndex,endIndex]
stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
end associate
@ -324,13 +332,15 @@ end subroutine phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine phenopowerlaw_dotState(Mp,ph,en)
module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal) :: &
xi_sl_sat_offset,&
@ -340,28 +350,32 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en)
right_SlipSlip
associate(prm => param(ph), stt => state(ph), dot => dotState(ph))
associate(prm => param(ph), stt => state(ph), &
dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), &
dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), &
dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), &
dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2)))
call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en))
dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
call kinetics_tw(Mp,ph,en,dot_gamma_tw)
sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, &
1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
dot%xi_sl(:,en) = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
* matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot%gamma_tw(:,en))
dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
* matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot_gamma_tw)
dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
* matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) &
+ prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en))
dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
* matmul(prm%h_tw_sl,dot_gamma_sl) &
+ prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw)
end associate
end subroutine phenopowerlaw_dotState
end function phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------

30
src/phase_thermal.f90
View File

@ -254,6 +254,36 @@ function integrateThermalState(Delta_t, ph,en) result(broken)
end function integrateThermalState
module subroutine thermal_restartWrite(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph
integer :: so
do so = 1,thermal_Nsources(ph)
call HDF5_write(thermalState(ph)%p(so)%state,groupHandle,'omega_thermal')
enddo
end subroutine thermal_restartWrite
module subroutine thermal_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph
integer :: so
do so = 1,thermal_Nsources(ph)
call HDF5_read(thermalState(ph)%p(so)%state0,groupHandle,'omega_thermal')
enddo
end subroutine thermal_restartRead
module subroutine thermal_forward()
integer :: ph, so

179
src/polynomials.f90 Normal file
View File

@ -0,0 +1,179 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Polynomial representation for variable data
!--------------------------------------------------------------------------------------------------
module polynomials
use prec
use IO
use YAML_parse
use YAML_types
implicit none
private
type, public :: tPolynomial
real(pReal), dimension(:), allocatable :: coef
real(pReal) :: x_ref
contains
procedure, public :: at => eval
procedure, public :: der1_at => eval_der1
end type tPolynomial
interface polynomial
module procedure polynomial_from_dict
module procedure polynomial_from_coef
end interface polynomial
public :: &
polynomial, &
polynomials_init
contains
!--------------------------------------------------------------------------------------------------
!> @brief Run self-test.
!--------------------------------------------------------------------------------------------------
subroutine polynomials_init()
print'(/,1x,a)', '<<<+- polynomials init -+>>>'; flush(IO_STDOUT)
call selfTest()
end subroutine polynomials_init
!--------------------------------------------------------------------------------------------------
!> @brief Initialize a Polynomial from Coefficients.
!--------------------------------------------------------------------------------------------------
function polynomial_from_coef(coef,x_ref) result(p)
real(pReal), dimension(:), intent(in) :: coef
real(pReal), intent(in) :: x_ref
type(tPolynomial) :: p
allocate(p%coef(0:size(coef)-1),source=coef) ! should be zero based
p%x_ref = x_ref
end function polynomial_from_coef
!--------------------------------------------------------------------------------------------------
!> @brief Initialize a Polynomial from a Dictionary with Coefficients.
!--------------------------------------------------------------------------------------------------
function polynomial_from_dict(dict,y,x) result(p)
type(tDict), intent(in) :: dict
character(len=*), intent(in) :: y, x
type(tPolynomial) :: p
real(pReal), dimension(:), allocatable :: coef
real(pReal) :: x_ref
allocate(coef(1),source=dict%get_asFloat(y))
if (dict%contains(y//','//x)) then
x_ref = dict%get_asFloat(x//'_ref')
coef = [coef,dict%get_asFloat(y//','//x)]
if (dict%contains(y//','//x//'^2')) then
coef = [coef,dict%get_asFloat(y//','//x//'^2')]
end if
else
x_ref = huge(0.0_pReal) ! Simplify debugging
end if
p = Polynomial(coef,x_ref)
end function polynomial_from_dict
!--------------------------------------------------------------------------------------------------
!> @brief Evaluate a Polynomial.
!--------------------------------------------------------------------------------------------------
pure function eval(self,x) result(y)
class(tPolynomial), intent(in) :: self
real(pReal), intent(in) :: x
real(pReal) :: y
integer :: i
y = self%coef(0)
do i = 1, ubound(self%coef,1)
y = y + self%coef(i) * (x-self%x_ref)**i
enddo
end function eval
!--------------------------------------------------------------------------------------------------
!> @brief Evaluate a first derivative of Polynomial.
!--------------------------------------------------------------------------------------------------
pure function eval_der1(self,x) result(y)
class(tPolynomial), intent(in) :: self
real(pReal), intent(in) :: x
real(pReal) :: y
integer :: i
y = 0.0_pReal
do i = 1, ubound(self%coef,1)
y = y + real(i,pReal)*self%coef(i) * (x-self%x_ref)**(i-1)
enddo
end function eval_der1
!--------------------------------------------------------------------------------------------------
!> @brief Check correctness of polynomical functionality.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
type(tPolynomial) :: p1, p2
real(pReal), dimension(3) :: coef
real(pReal) :: x_ref, x
class(tNode), pointer :: dict
character(len=pStringLen), dimension(3) :: coef_s
character(len=pStringLen) :: x_ref_s, x_s, YAML_s
call random_number(coef)
call random_number(x_ref)
call random_number(x)
coef = coef*10_pReal -0.5_pReal
x_ref = x_ref*10_pReal -0.5_pReal
x = x*10_pReal -0.5_pReal
p1 = polynomial(coef,x_ref)
if (dNeq(p1%at(x_ref),coef(1))) error stop 'polynomial: @ref'
write(coef_s(1),*) coef(1)
write(coef_s(2),*) coef(2)
write(coef_s(3),*) coef(3)
write(x_ref_s,*) x_ref
write(x_s,*) x
YAML_s = 'C: '//trim(adjustl(coef_s(1)))//IO_EOL//&
'C,T: '//trim(adjustl(coef_s(2)))//IO_EOL//&
'C,T^2: '//trim(adjustl(coef_s(3)))//IO_EOL//&
'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL
Dict => YAML_parse_str(trim(YAML_s))
p2 = polynomial(dict%asDict(),'C','T')
if (dNeq(p1%at(x),p2%at(x),1.0e-10_pReal)) error stop 'polynomials: init'
p1 = polynomial(coef*[0.0_pReal,1.0_pReal,0.0_pReal],x_ref)
if (dNeq(p1%at(x_ref+x),-p1%at(x_ref-x),1.0e-10_pReal)) error stop 'polynomials: eval(odd)'
if (dNeq(p1%der1_at(x),p1%der1_at(5.0_pReal*x),1.0e-10_pReal)) error stop 'polynomials: eval_der(odd)'
p1 = polynomial(coef*[0.0_pReal,0.0_pReal,1.0_pReal],x_ref)
if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval(even)'
if (dNeq(p1%der1_at(x_ref+x),-p1%der1_at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval_der(even)'
end subroutine selfTest
end module polynomials

2
src/prec.f90
View File

@ -45,6 +45,8 @@ module prec
state, & !< state
dotState, & !< rate of state change
deltaState !< increment of state change
real(pReal), pointer, dimension(:,:) :: &
deltaState2
end type
type, extends(tState) :: tPlasticState

98
src/rotations.f90
View File

@ -55,7 +55,7 @@ module rotations
real(pReal), parameter :: P = -1.0_pReal !< parameter for orientation conversion.
type, public :: rotation
type, public :: tRotation
real(pReal), dimension(4) :: q
contains
procedure, public :: asQuaternion
@ -78,10 +78,9 @@ module rotations
procedure, public :: rotStiffness
procedure, public :: misorientation
procedure, public :: standardize
end type rotation
end type tRotation
real(pReal), parameter :: &
SPI = sqrt(PI), &
PREF = sqrt(6.0_pReal/PI), &
A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), &
AP = PI**(2.0_pReal/3.0_pReal), &
@ -118,8 +117,8 @@ end subroutine rotations_init
!---------------------------------------------------------------------------------------------------
pure function asQuaternion(self)
class(rotation), intent(in) :: self
real(pReal), dimension(4) :: asQuaternion
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asQuaternion
asQuaternion = self%q
@ -127,8 +126,8 @@ end function asQuaternion
!---------------------------------------------------------------------------------------------------
pure function asEulers(self)
class(rotation), intent(in) :: self
real(pReal), dimension(3) :: asEulers
class(tRotation), intent(in) :: self
real(pReal), dimension(3) :: asEulers
asEulers = qu2eu(self%q)
@ -136,8 +135,8 @@ end function asEulers
!---------------------------------------------------------------------------------------------------
pure function asAxisAngle(self)
class(rotation), intent(in) :: self
real(pReal), dimension(4) :: asAxisAngle
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asAxisAngle
asAxisAngle = qu2ax(self%q)
@ -145,8 +144,8 @@ end function asAxisAngle
!---------------------------------------------------------------------------------------------------
pure function asMatrix(self)
class(rotation), intent(in) :: self
real(pReal), dimension(3,3) :: asMatrix
class(tRotation), intent(in) :: self
real(pReal), dimension(3,3) :: asMatrix
asMatrix = qu2om(self%q)
@ -154,8 +153,8 @@ end function asMatrix
!---------------------------------------------------------------------------------------------------
pure function asRodrigues(self)
class(rotation), intent(in) :: self
real(pReal), dimension(4) :: asRodrigues
class(tRotation), intent(in) :: self
real(pReal), dimension(4) :: asRodrigues
asRodrigues = qu2ro(self%q)
@ -163,8 +162,8 @@ end function asRodrigues
!---------------------------------------------------------------------------------------------------
pure function asHomochoric(self)
class(rotation), intent(in) :: self
real(pReal), dimension(3) :: asHomochoric
class(tRotation), intent(in) :: self
real(pReal), dimension(3) :: asHomochoric
asHomochoric = qu2ho(self%q)
@ -175,7 +174,7 @@ end function asHomochoric
!---------------------------------------------------------------------------------------------------
subroutine fromQuaternion(self,qu)
class(rotation), intent(out) :: self
class(tRotation), intent(out) :: self
real(pReal), dimension(4), intent(in) :: qu
if (dNeq(norm2(qu),1.0_pReal,1.0e-8_pReal)) call IO_error(402,ext_msg='fromQuaternion')
@ -186,7 +185,7 @@ end subroutine fromQuaternion
!---------------------------------------------------------------------------------------------------
subroutine fromEulers(self,eu,degrees)
class(rotation), intent(out) :: self
class(tRotation), intent(out) :: self
real(pReal), dimension(3), intent(in) :: eu
logical, intent(in), optional :: degrees
@ -198,7 +197,7 @@ subroutine fromEulers(self,eu,degrees)
Eulers = merge(eu*INRAD,eu,degrees)
endif
if (any(Eulers<0.0_pReal) .or. any(Eulers>2.0_pReal*PI) .or. Eulers(2) > PI) &
if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) &
call IO_error(402,ext_msg='fromEulers')
self%q = eu2qu(Eulers)
@ -207,7 +206,7 @@ end subroutine fromEulers
!---------------------------------------------------------------------------------------------------
subroutine fromAxisAngle(self,ax,degrees,P)
class(rotation), intent(out) :: self
class(tRotation), intent(out) :: self
real(pReal), dimension(4), intent(in) :: ax
logical, intent(in), optional :: degrees
integer, intent(in), optional :: P
@ -237,7 +236,7 @@ end subroutine fromAxisAngle
!---------------------------------------------------------------------------------------------------
subroutine fromMatrix(self,om)
class(rotation), intent(out) :: self
class(tRotation), intent(out) :: self
real(pReal), dimension(3,3), intent(in) :: om
if (dNeq(math_det33(om),1.0_pReal,tol=1.0e-5_pReal)) &
@ -254,10 +253,10 @@ end subroutine fromMatrix
!---------------------------------------------------------------------------------------------------
pure elemental function rotRot__(self,R) result(rRot)
type(rotation) :: rRot
class(rotation), intent(in) :: self,R
type(tRotation) :: rRot
class(tRotation), intent(in) :: self,R
rRot = rotation(multiply_quaternion(self%q,R%q))
rRot = tRotation(multiply_quaternion(self%q,R%q))
call rRot%standardize()
end function rotRot__
@ -268,7 +267,7 @@ end function rotRot__
!---------------------------------------------------------------------------------------------------
pure elemental subroutine standardize(self)
class(rotation), intent(inout) :: self
class(tRotation), intent(inout) :: self
if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q
@ -282,7 +281,7 @@ end subroutine standardize
pure function rotVector(self,v,active) result(vRot)
real(pReal), dimension(3) :: vRot
class(rotation), intent(in) :: self
class(tRotation), intent(in) :: self
real(pReal), intent(in), dimension(3) :: v
logical, intent(in), optional :: active
@ -318,7 +317,7 @@ end function rotVector
pure function rotTensor2(self,T,active) result(tRot)
real(pReal), dimension(3,3) :: tRot
class(rotation), intent(in) :: self
class(tRotation), intent(in) :: self
real(pReal), intent(in), dimension(3,3) :: T
logical, intent(in), optional :: active
@ -347,7 +346,7 @@ end function rotTensor2
pure function rotTensor4(self,T,active) result(tRot)
real(pReal), dimension(3,3,3,3) :: tRot
class(rotation), intent(in) :: self
class(tRotation), intent(in) :: self
real(pReal), intent(in), dimension(3,3,3,3) :: T
logical, intent(in), optional :: active
@ -379,7 +378,7 @@ end function rotTensor4
pure function rotStiffness(self,C,active) result(cRot)
real(pReal), dimension(6,6) :: cRot
class(rotation), intent(in) :: self
class(tRotation), intent(in) :: self
real(pReal), intent(in), dimension(6,6) :: C
logical, intent(in), optional :: active
@ -416,8 +415,8 @@ end function rotStiffness
!---------------------------------------------------------------------------------------------------
pure elemental function misorientation(self,other)
type(rotation) :: misorientation
class(rotation), intent(in) :: self, other
type(tRotation) :: misorientation
class(tRotation), intent(in) :: self, other
misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q))
@ -480,7 +479,7 @@ pure function qu2eu(qu) result(eu)
atan2( 2.0_pReal*chi, q03-q12 ), &
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
endif degenerated
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU])
end function qu2eu
@ -628,7 +627,7 @@ pure function om2eu(om) result(eu)
eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ]
end if
where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU])
end function om2eu
@ -1209,7 +1208,7 @@ pure function ho2cu(ho) result(cu)
else special
q2 = qxy + maxval(abs(xyz2))**2
sq2 = sqrt(q2)
q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
q = (BETA/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy
Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], &
[ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], &
@ -1217,7 +1216,7 @@ pure function ho2cu(ho) result(cu)
endif special
! inverse M_1
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / PREF ]/SC
! reverse the coordinates back to order according to the original pyramid number
cu = xyz1(p(:,2))
@ -1323,32 +1322,32 @@ pure function cu2ho(cu) result(ho)
else center
! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cu)
XYZ = cu(p(:,1)) * sc
XYZ = cu(p(:,1)) * SC
! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then
LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ]
LamXYZ = [ 0.0_pReal, 0.0_pReal, PREF * XYZ(3) ]
else special
order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ
q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger
c = cos(q)
s = sin(q)
q = prek * XYZ(order(2))/ sqrt(R2-c)
q = PREK * XYZ(order(2))/ sqrt(R2-c)
T = [ (R2*c - 1.0), R2 * s] * q
! transform to sphere grid (inverse Lambert)
! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
c = sum(T**2)
s = Pi * c/(24.0*XYZ(3)**2)
c = sPi * c / sqrt(24.0_pReal) / XYZ(3)
s = PI * c/(24.0*XYZ(3)**2)
c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3)
q = sqrt( 1.0 - s )
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ]
endif special
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ]
end if special
! reverse the coordinates back to order according to the original pyramid number
ho = LamXYZ(p(:,2))
endif center
end if center
end function cu2ho
@ -1416,7 +1415,7 @@ end function conjugate_quaternion
!--------------------------------------------------------------------------------------------------
subroutine selfTest()
type(rotation) :: R
type(tRotation) :: R
real(pReal), dimension(4) :: qu, ax, ro
real(pReal), dimension(3) :: x, eu, ho, v3
real(pReal), dimension(3,3) :: om, t33
@ -1437,7 +1436,7 @@ subroutine selfTest()
elseif(i==2) then
qu = eu2qu([0.0_pReal,0.0_pReal,0.0_pReal])
elseif(i==3) then
qu = eu2qu([2.0_pReal*PI,PI,2.0_pReal*PI])
qu = eu2qu([TAU,PI,TAU])
elseif(i==4) then
qu = [0.0_pReal,0.0_pReal,1.0_pReal,0.0_pReal]
elseif(i==5) then
@ -1448,10 +1447,10 @@ subroutine selfTest()
call random_number(x)
A = sqrt(x(3))
B = sqrt(1-0_pReal -x(3))
qu = [cos(2.0_pReal*PI*x(1))*A,&
sin(2.0_pReal*PI*x(2))*B,&
cos(2.0_pReal*PI*x(2))*B,&
sin(2.0_pReal*PI*x(1))*A]
qu = [cos(TAU*x(1))*A,&
sin(TAU*x(2))*B,&
cos(TAU*x(2))*B,&
sin(TAU*x(1))*A]
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif
@ -1504,7 +1503,8 @@ subroutine selfTest()
call random_number(C)
C = C+transpose(C)
if (any(dNeq(R%rotStiffness(C),math_3333toVoigt66(R%rotate(math_Voigt66to3333(C))),1.0e-12_pReal))) &
if (any(dNeq(R%rotStiffness(C), &
math_3333toVoigt66_stiffness(R%rotate(math_Voigt66to3333_stiffness(C))),1.0e-12_pReal))) &
error stop 'rotStiffness'
call R%fromQuaternion(qu * (1.0_pReal + merge(+5.e-9_pReal,-5.e-9_pReal, mod(i,2) == 0))) ! allow reasonable tolerance for ASCII/YAML

21
src/system_routines.f90
View File

@ -150,14 +150,14 @@ function getUserName()
getUserName = c_f_string(getUserName_Cstring)
else
getUserName = 'n/a (Error!)'
endif
end if
end function getUserName
!--------------------------------------------------------------------------------------------------
!> @brief convert C string to Fortran string
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string
!> @brief Convert C string to Fortran string.
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string.
!--------------------------------------------------------------------------------------------------
pure function c_f_string(c_string) result(f_string)
@ -174,28 +174,23 @@ pure function c_f_string(c_string) result(f_string)
else
f_string = f_string(:i-1)
exit
endif
enddo arrayToString
end if
end do arrayToString
end function c_f_string
!--------------------------------------------------------------------------------------------------
!> @brief convert Fortran string to C string
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string
!> @brief Convert Fortran string to C string.
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string.
!--------------------------------------------------------------------------------------------------
pure function f_c_string(f_string) result(c_string)
character(len=*), intent(in) :: f_string
character(kind=C_CHAR), dimension(len_trim(f_string)+1) :: c_string
integer :: i
do i=1,len_trim(f_string)
c_string(i)=f_string(i:i)
enddo
c_string(len_trim(f_string)+1) = C_NULL_CHAR
c_string = transfer(trim(f_string)//C_NULL_CHAR,c_string,size=size(c_string))
end function f_c_string