Merge remote-tracking branch 'origin/development' into MSC-Version

This commit is contained in:
Martin Diehl 2021-12-01 15:13:46 +01:00
commit 5abfe3c214
24 changed files with 440 additions and 239 deletions

@ -1 +1 @@
Subproject commit 76bb51348de75207d483d369628670e5ae51dca9
Subproject commit 02609955e53c0a5fb6cee5753419fb1ba1b9da2a

View File

@ -6,7 +6,7 @@ references:
output: [xi_sl, xi_tw]
N_sl: [3, 3, 0, 6, 0, 6] # basal, 1. prism, -, 1. pyr<a>, -, 2. pyr<c+a>
N_sl: [3, 3, 0, 6, 0, 6] # basal, prism, -, 1. pyr<a>, -, 2. pyr<c+a>
N_tw: [6, 0, 6] # tension, -, compression
xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6]
@ -23,7 +23,16 @@ f_sat_sl-tw: 10.0
h_0_sl-sl: 500.0e+6
h_0_tw-tw: 50.0e+6
h_0_tw-sl: 150.0e+6
h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0,
+1.0, 1.0] # unused entries are indicated by -1.0
h_tw-tw: [+1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0,
-1.0, 1.0] # unused entries are indicated by -1.0
h_tw-sl: [+1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+1.0, -1.0, 1.0, -1.0] # unused entries are indicated by -1.0
h_sl-tw: [+1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+1.0, -1.0, 1.0] # unused entries are indicated by -1.0

View File

@ -8,7 +8,7 @@ references:
https://doi.org/10.1016/j.actamat.2017.05.015
output: [gamma_sl]
N_sl: [3, 3, 0, 0, 12] # basal, 1. prism, -, -, 2. pyr<c+a>
N_sl: [3, 3, 0, 0, 12] # basal, 1. prism, -, -, 1. pyr<c+a>
n_sl: 20
a_sl: 2.0
dot_gamma_0_sl: 0.001
@ -20,4 +20,6 @@ xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6]
# L. Wang et al. :
# xi_0_sl: [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6]
h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0] # unused entries are indicated by -1.0

View File

@ -1 +1 @@
v3.0.0-alpha5-155-gbf76d9f3a
v3.0.0-alpha5-210-g7e7098baf

View File

@ -2,6 +2,8 @@ import os
import json
import functools
import colorsys
from pathlib import Path
from typing import Sequence, Union, TextIO
import numpy as np
import matplotlib as mpl
@ -14,9 +16,9 @@ from PIL import Image
from . import util
from . import Table
_eps = 216./24389.
_kappa = 24389./27.
_ref_white = np.array([.95047, 1.00000, 1.08883]) # Observer = 2, Illuminant = D65
_EPS = 216./24389.
_KAPPA = 24389./27.
_REF_WHITE = np.array([.95047, 1.00000, 1.08883]) # Observer = 2, Illuminant = D65
# ToDo (if needed)
# - support alpha channel (paraview/ASCII/input)
@ -39,20 +41,20 @@ class Colormap(mpl.colors.ListedColormap):
"""
def __add__(self,other):
def __add__(self, other: 'Colormap') -> 'Colormap':
"""Concatenate."""
return Colormap(np.vstack((self.colors,other.colors)),
f'{self.name}+{other.name}')
def __iadd__(self,other):
def __iadd__(self, other: 'Colormap') -> 'Colormap':
"""Concatenate (in-place)."""
return self.__add__(other)
def __invert__(self):
def __invert__(self) -> 'Colormap':
"""Reverse."""
return self.reversed()
def __repr__(self):
def __repr__(self) -> str:
"""Show as matplotlib figure."""
fig = plt.figure(self.name,figsize=(5,.5))
ax1 = fig.add_axes([0, 0, 1, 1])
@ -64,7 +66,11 @@ class Colormap(mpl.colors.ListedColormap):
@staticmethod
def from_range(low,high,name='DAMASK colormap',N=256,model='rgb'):
def from_range(low: Sequence[float],
high: Sequence[float],
name: str = 'DAMASK colormap',
N: int = 256,
model: str = 'rgb') -> 'Colormap':
"""
Create a perceptually uniform colormap between given (inclusive) bounds.
@ -145,7 +151,7 @@ class Colormap(mpl.colors.ListedColormap):
@staticmethod
def from_predefined(name,N=256):
def from_predefined(name: str, N: int = 256) -> 'Colormap':
"""
Select from a set of predefined colormaps.
@ -185,7 +191,10 @@ class Colormap(mpl.colors.ListedColormap):
return Colormap.from_range(definition['low'],definition['high'],name,N)
def shade(self,field,bounds=None,gap=None):
def shade(self,
field: np.ndarray,
bounds: Sequence[float] = None,
gap: float = None) -> Image:
"""
Generate PIL image of 2D field using colormap.
@ -226,7 +235,7 @@ class Colormap(mpl.colors.ListedColormap):
mode='RGBA')
def reversed(self,name=None):
def reversed(self, name: str = None) -> 'Colormap':
"""
Reverse.
@ -251,7 +260,7 @@ class Colormap(mpl.colors.ListedColormap):
return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name)
def _get_file_handle(self,fname,extension):
def _get_file_handle(self, fname: Union[TextIO, str, Path, None], suffix: str) -> TextIO:
"""
Provide file handle.
@ -259,8 +268,7 @@ class Colormap(mpl.colors.ListedColormap):
----------
fname : file, str, pathlib.Path, or None
Filename or filehandle, will be name of the colormap+extension if None.
extension: str
suffix: str
Extension of the filename.
Returns
@ -270,17 +278,14 @@ class Colormap(mpl.colors.ListedColormap):
"""
if fname is None:
fhandle = open(self.name.replace(' ','_')+'.'+extension,'w',newline='\n')
return open(self.name.replace(' ','_')+suffix, 'w', newline='\n')
elif isinstance(fname, (str, Path)):
return open(fname, 'w', newline='\n')
else:
try:
fhandle = open(fname,'w',newline='\n')
except TypeError:
fhandle = fname
return fhandle
return fname
def save_paraview(self,fname=None):
def save_paraview(self, fname: Union[TextIO, str, Path] = None):
"""
Save as JSON file for use in Paraview.
@ -303,12 +308,12 @@ class Colormap(mpl.colors.ListedColormap):
'RGBPoints':colors
}]
fhandle = self._get_file_handle(fname,'json')
fhandle = self._get_file_handle(fname,'.json')
json.dump(out,fhandle,indent=4)
fhandle.write('\n')
def save_ASCII(self,fname=None):
def save_ASCII(self, fname: Union[TextIO, str, Path] = None):
"""
Save as ASCII file.
@ -321,10 +326,10 @@ class Colormap(mpl.colors.ListedColormap):
"""
labels = {'RGBA':4} if self.colors.shape[1] == 4 else {'RGB': 3}
t = Table(self.colors,labels,f'Creator: {util.execution_stamp("Colormap")}')
t.save(self._get_file_handle(fname,'txt'))
t.save(self._get_file_handle(fname,'.txt'))
def save_GOM(self,fname=None):
def save_GOM(self, fname: Union[TextIO, str, Path] = None):
"""
Save as ASCII file for use in GOM Aramis.
@ -342,10 +347,10 @@ class Colormap(mpl.colors.ListedColormap):
+ ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(int))]) \
+ '\n'
self._get_file_handle(fname,'legend').write(GOM_str)
self._get_file_handle(fname,'.legend').write(GOM_str)
def save_gmsh(self,fname=None):
def save_gmsh(self, fname: Union[TextIO, str, Path] = None):
"""
Save as ASCII file for use in gmsh.
@ -360,11 +365,13 @@ class Colormap(mpl.colors.ListedColormap):
gmsh_str = 'View.ColorTable = {\n' \
+'\n'.join([f'{c[0]},{c[1]},{c[2]},' for c in self.colors[:,:3]*255]) \
+'\n}\n'
self._get_file_handle(fname,'msh').write(gmsh_str)
self._get_file_handle(fname,'.msh').write(gmsh_str)
@staticmethod
def _interpolate_msh(frac,low,high):
def _interpolate_msh(frac,
low: np.ndarray,
high: np.ndarray) -> np.ndarray:
"""
Interpolate in Msh color space.
@ -441,31 +448,31 @@ class Colormap(mpl.colors.ListedColormap):
@staticmethod
def _hsv2rgb(hsv):
def _hsv2rgb(hsv: np.ndarray) -> np.ndarray:
"""H(ue) S(aturation) V(alue) to R(red) G(reen) B(lue)."""
return np.array(colorsys.hsv_to_rgb(hsv[0]/360.,hsv[1],hsv[2]))
@staticmethod
def _rgb2hsv(rgb):
def _rgb2hsv(rgb: np.ndarray) -> np.ndarray:
"""R(ed) G(reen) B(lue) to H(ue) S(aturation) V(alue)."""
h,s,v = colorsys.rgb_to_hsv(rgb[0],rgb[1],rgb[2])
return np.array([h*360,s,v])
@staticmethod
def _hsl2rgb(hsl):
def _hsl2rgb(hsl: np.ndarray) -> np.ndarray:
"""H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue)."""
return np.array(colorsys.hls_to_rgb(hsl[0]/360.,hsl[2],hsl[1]))
@staticmethod
def _rgb2hsl(rgb):
def _rgb2hsl(rgb: np.ndarray) -> np.ndarray:
"""R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance)."""
h,l,s = colorsys.rgb_to_hls(rgb[0],rgb[1],rgb[2])
return np.array([h*360,s,l])
@staticmethod
def _xyz2rgb(xyz):
def _xyz2rgb(xyz: np.ndarray) -> np.ndarray:
"""
CIE Xyz to R(ed) G(reen) B(lue).
@ -485,7 +492,7 @@ class Colormap(mpl.colors.ListedColormap):
return np.clip(rgb,0.,1.)
@staticmethod
def _rgb2xyz(rgb):
def _rgb2xyz(rgb: np.ndarray) -> np.ndarray:
"""
R(ed) G(reen) B(lue) to CIE Xyz.
@ -503,7 +510,7 @@ class Colormap(mpl.colors.ListedColormap):
@staticmethod
def _lab2xyz(lab,ref_white=None):
def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
"""
CIE Lab to CIE Xyz.
@ -516,13 +523,13 @@ class Colormap(mpl.colors.ListedColormap):
f_z = (lab[0]+16.)/116. - lab[2]/200.
return np.array([
f_x**3. if f_x**3. > _eps else (116.*f_x-16.)/_kappa,
((lab[0]+16.)/116.)**3 if lab[0]>_kappa*_eps else lab[0]/_kappa,
f_z**3. if f_z**3. > _eps else (116.*f_z-16.)/_kappa
])*(ref_white if ref_white is not None else _ref_white)
f_x**3. if f_x**3. > _EPS else (116.*f_x-16.)/_KAPPA,
((lab[0]+16.)/116.)**3 if lab[0]>_KAPPA*_EPS else lab[0]/_KAPPA,
f_z**3. if f_z**3. > _EPS else (116.*f_z-16.)/_KAPPA
])*(ref_white if ref_white is not None else _REF_WHITE)
@staticmethod
def _xyz2lab(xyz,ref_white=None):
def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
"""
CIE Xyz to CIE Lab.
@ -531,8 +538,8 @@ class Colormap(mpl.colors.ListedColormap):
http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
"""
ref_white = ref_white if ref_white is not None else _ref_white
f = np.where(xyz/ref_white > _eps,(xyz/ref_white)**(1./3.),(_kappa*xyz/ref_white+16.)/116.)
ref_white = ref_white if ref_white is not None else _REF_WHITE
f = np.where(xyz/ref_white > _EPS,(xyz/ref_white)**(1./3.),(_KAPPA*xyz/ref_white+16.)/116.)
return np.array([
116.0 * f[1] - 16.0,
@ -542,7 +549,7 @@ class Colormap(mpl.colors.ListedColormap):
@staticmethod
def _lab2msh(lab):
def _lab2msh(lab: np.ndarray) -> np.ndarray:
"""
CIE Lab to Msh.
@ -560,7 +567,7 @@ class Colormap(mpl.colors.ListedColormap):
])
@staticmethod
def _msh2lab(msh):
def _msh2lab(msh: np.ndarray) -> np.ndarray:
"""
Msh to CIE Lab.
@ -577,29 +584,29 @@ class Colormap(mpl.colors.ListedColormap):
])
@staticmethod
def _lab2rgb(lab):
def _lab2rgb(lab: np.ndarray) -> np.ndarray:
return Colormap._xyz2rgb(Colormap._lab2xyz(lab))
@staticmethod
def _rgb2lab(rgb):
def _rgb2lab(rgb: np.ndarray) -> np.ndarray:
return Colormap._xyz2lab(Colormap._rgb2xyz(rgb))
@staticmethod
def _msh2rgb(msh):
def _msh2rgb(msh: np.ndarray) -> np.ndarray:
return Colormap._lab2rgb(Colormap._msh2lab(msh))
@staticmethod
def _rgb2msh(rgb):
def _rgb2msh(rgb: np.ndarray) -> np.ndarray:
return Colormap._lab2msh(Colormap._rgb2lab(rgb))
@staticmethod
def _hsv2msh(hsv):
def _hsv2msh(hsv: np.ndarray) -> np.ndarray:
return Colormap._rgb2msh(Colormap._hsv2rgb(hsv))
@staticmethod
def _hsl2msh(hsl):
def _hsl2msh(hsl: np.ndarray) -> np.ndarray:
return Colormap._rgb2msh(Colormap._hsl2rgb(hsl))
@staticmethod
def _xyz2msh(xyz):
def _xyz2msh(xyz: np.ndarray) -> np.ndarray:
return Colormap._lab2msh(Colormap._xyz2lab(xyz))

View File

@ -77,6 +77,12 @@ class TestColormap:
# xyz2msh
assert np.allclose(Colormap._xyz2msh(xyz),msh,atol=1.e-6,rtol=0)
@pytest.mark.parametrize('low,high',[((0,0,0),(1,1,1)),
([0,0,0],[1,1,1])])
def test_from_range_types(self,low,high):
a = Colormap.from_range(low,high)
b = Colormap.from_range(np.array(low),np.array(high))
assert np.all(a.colors == b.colors)
@pytest.mark.parametrize('format',['ASCII','paraview','GOM','gmsh'])
@pytest.mark.parametrize('model',['rgb','hsv','hsl','xyz','lab','msh'])

View File

@ -390,7 +390,7 @@ class TestResult:
with open((ref_path/'export_VTK'/request.node.name).with_suffix('.md5'),'w') as f:
f.write(cur+'\n')
with open((ref_path/'export_VTK'/request.node.name).with_suffix('.md5')) as f:
assert cur == f.read()[:-1]
assert cur == f.read().strip('\n')
@pytest.mark.parametrize('mode',['point','cell'])
@pytest.mark.parametrize('output',[False,True])

View File

@ -4,6 +4,7 @@
!> @details List of files needed by MSC.Marc
!--------------------------------------------------------------------------------------------------
#include "parallelization.f90"
#include "constants.f90"
#include "IO.f90"
#include "YAML_types.f90"
#include "YAML_parse.f90"

15
src/constants.f90 Normal file
View File

@ -0,0 +1,15 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief physical constants
!--------------------------------------------------------------------------------------------------
module constants
use prec
implicit none
public
real(pReal), parameter :: &
T_ROOM = 300.0_pReal, & !< Room temperature in K
K_B = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
end module constants

View File

@ -223,7 +223,11 @@ program DAMASK_grid
loadCases(l)%r = step_discretization%get_asFloat('r', defaultVal= 1.0_pReal)
loadCases(l)%f_restart = load_step%get_asInt('f_restart', defaultVal=huge(0))
if (load_step%get_asString('f_out',defaultVal='n/a') == 'none') then
loadCases(l)%f_out = huge(0)
else
loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1)
end if
loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1)
reportAndCheck: if (worldrank == 0) then
@ -276,6 +280,7 @@ program DAMASK_grid
endif
print'(2x,a,1x,f0.3)', 't:', loadCases(l)%t
print'(2x,a,1x,i0)', 'N:', loadCases(l)%N
if (loadCases(l)%f_out < huge(0)) &
print'(2x,a,1x,i0)', 'f_out:', loadCases(l)%f_out
if (loadCases(l)%f_restart < huge(0)) &
print'(2x,a,1x,i0)', 'f_restart:', loadCases(l)%f_restart

View File

@ -697,9 +697,9 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! applying boundary conditions
diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + &
C_volAvg(2,2,2,2)/delta(2)**2.0_pReal + &
C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ
diag = (C_volAvg(1,1,1,1)/delta(1)**2 + &
C_volAvg(2,2,2,2)/delta(2)**2 + &
C_volAvg(3,3,3,3)/delta(3)**2)*detJ
call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr)
CHKERRQ(ierr)
call DMGetGlobalVector(da_local,coordinates,ierr); CHKERRQ(ierr)

View File

@ -578,20 +578,20 @@ real(pReal) function utilities_divergenceRMS()
do k = 1, grid3; do j = 1, grid(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. do not take square root and square again
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal)& ! --> sum squared L_2 norm of vector
+ 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
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2) & ! --> sum squared L_2 norm of vector
+sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),&
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal))
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)
+ sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) &
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) &
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
+ sum( real(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal) &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2) &
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal)
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
call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
@ -630,7 +630,7 @@ real(pReal) function utilities_curlRMS()
-tensorField_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2))
enddo
utilities_curlRMS = utilities_curlRMS &
+2.0_pReal*sum(curl_fourier%re**2.0_pReal+curl_fourier%im**2.0_pReal) ! Has somewhere a conj. complex counterpart. Therefore count it twice.
+2.0_pReal*sum(curl_fourier%re**2+curl_fourier%im**2) ! Has somewhere a conj. complex counterpart. Therefore count it twice.
enddo
do l = 1, 3
curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2) &
@ -641,7 +641,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.0_pReal + curl_fourier%im**2.0_pReal) ! 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 grid(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))
@ -651,7 +651,7 @@ 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.0_pReal + curl_fourier%im**2.0_pReal) ! 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 grid(1) /= 1)
enddo; enddo
call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
@ -836,13 +836,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
dPdF_min = huge(1.0_pReal)
dPdF_norm_min = huge(1.0_pReal)
do i = 1, product(grid(1:2))*grid3
if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) then
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.0_pReal)
dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)
endif
if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) then
if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)
dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)
endif
enddo

View File

@ -546,7 +546,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
do k = 1,3; do l = 1,3
nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient
end do; end do
nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor
nDefNorm = nDefNorm + nDef(i,j)**2 ! compute the norm of the mismatch tensor
end do; end do
nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity)
nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces)

View File

@ -219,18 +219,18 @@ module lattice
2, -1, -1, 0, 0, 1, -1, 0, &
-1, 2, -1, 0, -1, 0, 1, 0, &
-1, -1, 2, 0, 1, -1, 0, 0, &
! <-11.0>{11.0}/2nd order prismatic compound systems (plane normal independent of c/a-ratio)
! <-11.0>{11.0}/2. order prismatic compound systems (plane normal independent of c/a-ratio)
-1, 1, 0, 0, 1, 1, -2, 0, &
0, -1, 1, 0, -2, 1, 1, 0, &
1, 0, -1, 0, 1, -2, 1, 0, &
! <-1-1.0>{-11.1}/1st order pyramidal <a> systems (direction independent of c/a-ratio)
! <-1-1.0>{-11.1}/1. order pyramidal <a> systems (direction independent of c/a-ratio)
-1, 2, -1, 0, 1, 0, -1, 1, &
-2, 1, 1, 0, 0, 1, -1, 1, &
-1, -1, 2, 0, -1, 1, 0, 1, &
1, -2, 1, 0, -1, 0, 1, 1, &
2, -1, -1, 0, 0, -1, 1, 1, &
1, 1, -2, 0, 1, -1, 0, 1, &
! <11.3>{-10.1}/1st order pyramidal <c+a> systems (direction independent of c/a-ratio)
! <11.3>{-10.1}/1. order pyramidal <c+a> systems (direction independent of c/a-ratio)
-2, 1, 1, 3, 1, 0, -1, 1, &
-1, -1, 2, 3, 1, 0, -1, 1, &
-1, -1, 2, 3, 0, 1, -1, 1, &
@ -243,7 +243,7 @@ module lattice
-1, 2, -1, 3, 0, -1, 1, 1, &
-1, 2, -1, 3, 1, -1, 0, 1, &
-2, 1, 1, 3, 1, -1, 0, 1, &
! <11.3>{-1-1.2}/2nd order pyramidal <c+a> systems
! <11.3>{-1-1.2}/2. order pyramidal <c+a> systems
-1, -1, 2, 3, 1, 1, -2, 2, &
1, -2, 1, 3, -1, 2, -1, 2, &
2, -1, -1, 3, -2, 1, 1, 2, &
@ -405,7 +405,7 @@ module lattice
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @brief Module initialization
!--------------------------------------------------------------------------------------------------
subroutine lattice_init
@ -417,7 +417,7 @@ end subroutine lattice_init
!--------------------------------------------------------------------------------------------------
!> @brief characteristic shear for twinning
!> @brief Characteristic shear for twinning
!--------------------------------------------------------------------------------------------------
function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear)
@ -473,13 +473,13 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
p = sum(HEX_NTWINSYSTEM(1:f-1))+s
select case(HEX_SHEARTWIN(p)) ! from Christian & Mahajan 1995 p.29
case (1) ! <-10.1>{10.2}
characteristicShear(a) = (3.0_pReal-cOverA**2.0_pReal)/sqrt(3.0_pReal)/CoverA
characteristicShear(a) = (3.0_pReal-cOverA**2)/sqrt(3.0_pReal)/CoverA
case (2) ! <11.6>{-1-1.1}
characteristicShear(a) = 1.0_pReal/cOverA
case (3) ! <10.-2>{10.1}
characteristicShear(a) = (4.0_pReal*cOverA**2.0_pReal-9.0_pReal)/sqrt(48.0_pReal)/cOverA
characteristicShear(a) = (4.0_pReal*cOverA**2-9.0_pReal)/sqrt(48.0_pReal)/cOverA
case (4) ! <11.-3>{11.2}
characteristicShear(a) = 2.0_pReal*(cOverA**2.0_pReal-2.0_pReal)/3.0_pReal/cOverA
characteristicShear(a) = 2.0_pReal*(cOverA**2-2.0_pReal)/3.0_pReal/cOverA
end select
case default
call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice))
@ -491,7 +491,7 @@ end function lattice_characteristicShear_Twin
!--------------------------------------------------------------------------------------------------
!> @brief rotated elasticity matrices for twinning in 6x6-matrix notation
!> @brief Rotated elasticity matrices for twinning in 6x6-matrix notation
!--------------------------------------------------------------------------------------------------
function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
@ -522,14 +522,14 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
do i = 1, sum(Ntwin)
call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg?
lattice_C66_twin(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C66)))
lattice_C66_twin(1:6,1:6,i) = R%rotStiffness(C66)
enddo
end function lattice_C66_twin
!--------------------------------------------------------------------------------------------------
!> @brief rotated elasticity matrices for transformation in 6x6-matrix notation
!> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation
!--------------------------------------------------------------------------------------------------
function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
cOverA_trans,a_bcc,a_fcc)
@ -548,6 +548,8 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
!--------------------------------------------------------------------------------------------------
! elasticity matrix of the target phase in cube orientation
if (lattice_target == 'hP') then
! https://doi.org/10.1063/1.1663858 eq. (16), eq. (18), eq. (19)
! https://doi.org/10.1016/j.actamat.2016.07.032 eq. (47), eq. (48)
if (cOverA_trans < 1.0_pReal .or. cOverA_trans > 2.0_pReal) &
call IO_error(131,ext_msg='lattice_C66_trans: '//trim(lattice_target))
C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pReal*C_parent66(4,4))/2.0_pReal
@ -558,11 +560,11 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
C_bar66(1,4) = (C_parent66(1,1) - C_parent66(1,2) - 2.0_pReal*C_parent66(4,4)) /(3.0_pReal*sqrt(2.0_pReal))
C_target_unrotated66 = 0.0_pReal
C_target_unrotated66(1,1) = C_bar66(1,1) - C_bar66(1,4)**2.0_pReal/C_bar66(4,4)
C_target_unrotated66(1,2) = C_bar66(1,2) + C_bar66(1,4)**2.0_pReal/C_bar66(4,4)
C_target_unrotated66(1,1) = C_bar66(1,1) - C_bar66(1,4)**2/C_bar66(4,4)
C_target_unrotated66(1,2) = C_bar66(1,2) + C_bar66(1,4)**2/C_bar66(4,4)
C_target_unrotated66(1,3) = C_bar66(1,3)
C_target_unrotated66(3,3) = C_bar66(3,3)
C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2)))
C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2)))
C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP')
elseif (lattice_target == 'cI') then
if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) &
@ -581,14 +583,14 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
do i = 1,sum(Ntrans)
call R%fromMatrix(Q(1:3,1:3,i))
lattice_C66_trans(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C_target_unrotated66)))
lattice_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66)
enddo
end function lattice_C66_trans
!--------------------------------------------------------------------------------------------------
!> @brief non-Schmid projections for bcc with up to 6 coefficients
!> @brief Non-schmid projections for bcc with up to 6 coefficients
! Koester et al. 2012, Acta Materialia 60 (2012) 38943901, eq. (17)
! Gröger et al. 2008, Acta Materialia 56 (2008) 54125425, table 1
!--------------------------------------------------------------------------------------------------
@ -635,7 +637,7 @@ end function lattice_nonSchmidMatrix
!--------------------------------------------------------------------------------------------------
!> @brief slip-slip interaction matrix
!> @brief Slip-slip interaction matrix
!> details only active slip systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
@ -751,22 +753,23 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: &
HEX_INTERACTIONSLIPSLIP = reshape( [&
! basal prism 2. prism 1. pyr<a> 1. pyr<c+a> 2. pyr<c+a>
1, 2, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! -----> acting (forest)
2, 1, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! |
2, 1, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! | basal
2, 2, 1, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! |
! v
6, 6, 6, 4, 5, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & ! reacting (primary)
6, 6, 6, 5, 4, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, &
6, 6, 6, 5, 4, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & ! prism
6, 6, 6, 5, 5, 4, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, &
12,12,12, 11,11,11, 9,10,10, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, &
12,12,12, 11,11,11, 10, 9,10, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, &
12,12,12, 11,11,11, 10, 9,10, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, & ! 2. prism
12,12,12, 11,11,11, 10,10, 9, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, &
20,20,20, 19,19,19, 18,18,18, 16,17,17,17,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, &
20,20,20, 19,19,19, 18,18,18, 17,16,17,17,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, &
20,20,20, 19,19,19, 18,18,18, 17,17,16,17,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, &
20,20,20, 19,19,19, 18,18,18, 17,17,17,16,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, &
20,20,20, 19,19,19, 18,18,18, 17,17,17,16,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & ! 1. pyr<a>
20,20,20, 19,19,19, 18,18,18, 17,17,17,17,16,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, &
20,20,20, 19,19,19, 18,18,18, 17,17,17,17,17,16, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, &
@ -776,7 +779,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,25,26,26,26,26,26,26,26,26, 35,35,35,35,35,35, &
30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,25,26,26,26,26,26,26,26, 35,35,35,35,35,35, &
30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,25,26,26,26,26,26,26, 35,35,35,35,35,35, &
30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,25,26,26,26,26,26, 35,35,35,35,35,35, &
30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,25,26,26,26,26,26, 35,35,35,35,35,35, & ! 1. pyr<c+a>
30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,25,26,26,26,26, 35,35,35,35,35,35, &
30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,26,25,26,26,26, 35,35,35,35,35,35, &
30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,26,26,25,26,26, 35,35,35,35,35,35, &
@ -786,7 +789,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 36,37,37,37,37,37, &
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,36,37,37,37,37, &
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,36,37,37,37, &
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,36,37,37, &
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,36,37,37, & ! 2. pyr<c+a>
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,36,37, &
42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,37,36 &
],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for hex (onion peel naming scheme)
@ -883,7 +886,7 @@ end function lattice_interaction_SlipBySlip
!--------------------------------------------------------------------------------------------------
!> @brief twin-twin interaction matrix
!> @brief Twin-twin interaction matrix
!> details only active twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix)
@ -932,31 +935,32 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(
!< 3: other interaction
integer, dimension(HEX_NTWIN,HEX_NTWIN), parameter :: &
HEX_INTERACTIONTWINTWIN = reshape( [&
! <-10.1>{10.2} <11.6>{-1-1.1} <10.-2>{10.1} <11.-3>{11.2}
1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! -----> acting
2, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! |
2, 2, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! |
2, 2, 2, 1, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! v
2, 2, 2, 1, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! v <-10.1>{10.2}
2, 2, 2, 2, 1, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! reacting
2, 2, 2, 2, 2, 1, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, &
6, 6, 6, 6, 6, 6, 4, 5, 5, 5, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, &
6, 6, 6, 6, 6, 6, 5, 4, 5, 5, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, &
6, 6, 6, 6, 6, 6, 5, 5, 4, 5, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, &
6, 6, 6, 6, 6, 6, 5, 5, 5, 4, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, &
6, 6, 6, 6, 6, 6, 5, 5, 5, 4, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, & ! <11.6>{-1-1.1}
6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 4, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, &
6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 4, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, &
12,12,12,12,12,12, 11,11,11,11,11,11, 9,10,10,10,10,10, 15,15,15,15,15,15, &
12,12,12,12,12,12, 11,11,11,11,11,11, 10, 9,10,10,10,10, 15,15,15,15,15,15, &
12,12,12,12,12,12, 11,11,11,11,11,11, 10,10, 9,10,10,10, 15,15,15,15,15,15, &
12,12,12,12,12,12, 11,11,11,11,11,11, 10,10,10, 9,10,10, 15,15,15,15,15,15, &
12,12,12,12,12,12, 11,11,11,11,11,11, 10,10,10, 9,10,10, 15,15,15,15,15,15, & ! <10.-2>{10.1}
12,12,12,12,12,12, 11,11,11,11,11,11, 10,10,10,10, 9,10, 15,15,15,15,15,15, &
12,12,12,12,12,12, 11,11,11,11,11,11, 10,10,10,10,10, 9, 15,15,15,15,15,15, &
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 16,17,17,17,17,17, &
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,16,17,17,17,17, &
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,16,17,17,17, &
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17, &
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17, & ! <11.-3>{11.2}
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,16,17, &
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 &
],shape(HEX_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for hex
@ -981,7 +985,7 @@ end function lattice_interaction_TwinByTwin
!--------------------------------------------------------------------------------------------------
!> @brief trans-trans interaction matrix
!> @brief Trans-trans interaction matrix
!> details only active trans systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix)
@ -1023,7 +1027,7 @@ end function lattice_interaction_TransByTrans
!--------------------------------------------------------------------------------------------------
!> @brief slip-twin interaction matrix
!> @brief Slip-twin interaction matrix
!> details only active slip and twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix)
@ -1123,22 +1127,23 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
integer, dimension(HEX_NTWIN,HEX_NSLIP), parameter :: &
HEX_INTERACTIONSLIPTWIN = reshape( [&
! <-10.1>{10.2} <11.6>{-1-1.1} <10.-2>{10.1} <11.-3>{11.2}
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! ----> twin (acting)
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! |
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | basal
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! |
! v
5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & ! slip (reacting)
5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, &
5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & ! prism
5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, &
9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, &
9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, &
9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & ! 2.prism
9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, &
13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, &
13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, &
13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, &
13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, &
13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & ! 1. pyr<a>
13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, &
13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, &
@ -1148,7 +1153,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, &
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, &
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, &
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, &
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & ! 1. pyr<c+a>
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, &
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, &
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, &
@ -1158,7 +1163,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, &
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, &
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, &
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, &
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & ! 2. pyr<c+a>
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, &
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 &
],shape(HEX_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hex
@ -1186,7 +1191,7 @@ end function lattice_interaction_SlipByTwin
!--------------------------------------------------------------------------------------------------
!> @brief slip-trans interaction matrix
!> @brief Slip-trans interaction matrix
!> details only active slip and trans systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix)
@ -1239,7 +1244,7 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice)
!--------------------------------------------------------------------------------------------------
!> @brief twin-slip interaction matrix
!> @brief Twin-slip interaction matrix
!> details only active twin and slip systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix)
@ -1262,31 +1267,32 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) r
integer, dimension(HEX_NSLIP,HEX_NTWIN), parameter :: &
HEX_INTERACTIONTWINSLIP = reshape( [&
! basal prism 2. prism 1. pyr<a> 1. pyr<c+a> 2. pyr<c+a>
1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! ----> slip (acting)
1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! |
1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! |
1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! v
1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! v <-10.1>{10.2}
1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! twin (reacting)
1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, &
2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, &
2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, &
2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, &
2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, &
2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & ! <11.6>{-1-1.1}
2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, &
2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, &
3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, &
3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, &
3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, &
3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, &
3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & ! <10.-2>{10.1}
3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, &
3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, &
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, &
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, &
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, &
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, &
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & ! <11.-3>{11.2}
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, &
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 &
],shape(HEX_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hex
@ -1483,7 +1489,7 @@ end function lattice_SchmidMatrix_cleavage
!--------------------------------------------------------------------------------------------------
!> @brief slip direction of slip systems (|| b)
!> @brief Slip direction of slip systems (|| b)
!--------------------------------------------------------------------------------------------------
function lattice_slip_direction(Nslip,lattice,cOverA) result(d)
@ -1501,7 +1507,7 @@ end function lattice_slip_direction
!--------------------------------------------------------------------------------------------------
!> @brief normal direction of slip systems (|| n)
!> @brief Normal direction of slip systems (|| n)
!--------------------------------------------------------------------------------------------------
function lattice_slip_normal(Nslip,lattice,cOverA) result(n)
@ -1519,7 +1525,7 @@ end function lattice_slip_normal
!--------------------------------------------------------------------------------------------------
!> @brief transverse direction of slip systems (|| t = b x n)
!> @brief Transverse direction of slip systems (|| t = b x n)
!--------------------------------------------------------------------------------------------------
function lattice_slip_transverse(Nslip,lattice,cOverA) result(t)
@ -1537,7 +1543,7 @@ end function lattice_slip_transverse
!--------------------------------------------------------------------------------------------------
!> @brief labels of slip systems
!> @brief Labels of slip systems
!> details only active slip systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_labels_slip(Nslip,lattice) result(labels)
@ -1578,7 +1584,7 @@ end function lattice_labels_slip
!--------------------------------------------------------------------------------------------------
!> @brief return 3x3 tensor with symmetry according to given Bravais lattice
!> @brief Return 3x3 tensor with symmetry according to given Bravais lattice
!--------------------------------------------------------------------------------------------------
pure function lattice_symmetrize_33(T,lattice) result(T_sym)
@ -1605,7 +1611,7 @@ end function lattice_symmetrize_33
!--------------------------------------------------------------------------------------------------
!> @brief return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice
!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice
!> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962
!--------------------------------------------------------------------------------------------------
pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym)
@ -1651,7 +1657,7 @@ end function lattice_symmetrize_C66
!--------------------------------------------------------------------------------------------------
!> @brief labels of twin systems
!> @brief Labels for twin systems
!> details only active twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_labels_twin(Ntwin,lattice) result(labels)
@ -1689,7 +1695,7 @@ end function lattice_labels_twin
!--------------------------------------------------------------------------------------------------
!> @brief projection of the transverse direction onto the slip plane
!> @brief Projection of the transverse direction onto the slip plane
!> @details: This projection is used to calculate forest hardening for edge dislocations
!--------------------------------------------------------------------------------------------------
function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
@ -1713,7 +1719,7 @@ end function slipProjection_transverse
!--------------------------------------------------------------------------------------------------
!> @brief projection of the slip direction onto the slip plane
!> @brief Projection of the slip direction onto the slip plane
!> @details: This projection is used to calculate forest hardening for screw dislocations
!--------------------------------------------------------------------------------------------------
function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
@ -1779,7 +1785,7 @@ end function coordinateSystem_slip
!--------------------------------------------------------------------------------------------------
!> @brief populate reduced interaction matrix
!> @brief Populate reduced interaction matrix
!--------------------------------------------------------------------------------------------------
function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,values,matrix)
@ -1822,7 +1828,7 @@ end function buildInteraction
!--------------------------------------------------------------------------------------------------
!> @brief build a local coordinate system on slip, twin, trans, cleavage systems
!> @brief Build a local coordinate system on slip, twin, trans, cleavage systems
!> @details Order: Direction, plane (normal), and common perpendicular
!--------------------------------------------------------------------------------------------------
function buildCoordinateSystem(active,potential,system,lattice,cOverA)
@ -1889,7 +1895,7 @@ end function buildCoordinateSystem
!--------------------------------------------------------------------------------------------------
!> @brief helper function to define transformation systems
!> @brief Helper function to define transformation systems
! Needed to calculate Schmid matrix and rotated stiffness matrices.
! @details: set c/a = 0.0 for fcc -> bcc transformation
! set a_Xcc = 0.0 for fcc -> hex transformation
@ -2073,7 +2079,7 @@ end function getlabels
!--------------------------------------------------------------------------------------------------
!> @brief equivalent Poisson's ratio (ν)
!> @brief Equivalent Poisson's ratio (ν)
!> @details https://doi.org/10.1143/JPSJ.20.635
!--------------------------------------------------------------------------------------------------
function lattice_equivalent_nu(C,assumption) result(nu)
@ -2106,7 +2112,7 @@ end function lattice_equivalent_nu
!--------------------------------------------------------------------------------------------------
!> @brief equivalent shear modulus (μ)
!> @brief Equivalent shear modulus (μ)
!> @details https://doi.org/10.1143/JPSJ.20.635
!--------------------------------------------------------------------------------------------------
function lattice_equivalent_mu(C,assumption) result(mu)
@ -2135,7 +2141,7 @@ end function lattice_equivalent_mu
!--------------------------------------------------------------------------------------------------
!> @brief check correctness of some lattice functions
!> @brief Check correctness of some lattice functions.
!--------------------------------------------------------------------------------------------------
subroutine selfTest

View File

@ -854,7 +854,70 @@ end function math_66toSym3333
!--------------------------------------------------------------------------------------------------
!> @brief convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix
!> @brief Convert 6 Voigt stress vector into symmetric 3x3 tensor.
!--------------------------------------------------------------------------------------------------
pure function math_Voigt6to33_stress(sigma_tilde) result(sigma)
real(pReal), dimension(3,3) :: sigma
real(pReal), dimension(6), intent(in) :: sigma_tilde
sigma = reshape([sigma_tilde(1), sigma_tilde(6), sigma_tilde(5), &
sigma_tilde(6), sigma_tilde(2), sigma_tilde(4), &
sigma_tilde(5), sigma_tilde(4), sigma_tilde(3)],[3,3])
end function math_Voigt6to33_stress
!--------------------------------------------------------------------------------------------------
!> @brief Convert 6 Voigt strain vector into symmetric 3x3 tensor.
!--------------------------------------------------------------------------------------------------
pure function math_Voigt6to33_strain(epsilon_tilde) result(epsilon)
real(pReal), dimension(3,3) :: epsilon
real(pReal), dimension(6), intent(in) :: epsilon_tilde
epsilon = reshape([ epsilon_tilde(1), 0.5_pReal*epsilon_tilde(6), 0.5_pReal*epsilon_tilde(5), &
0.5_pReal*epsilon_tilde(6), epsilon_tilde(2), 0.5_pReal*epsilon_tilde(4), &
0.5_pReal*epsilon_tilde(5), 0.5_pReal*epsilon_tilde(4), epsilon_tilde(3)],[3,3])
end function math_Voigt6to33_strain
!--------------------------------------------------------------------------------------------------
!> @brief Convert 3x3 tensor into 6 Voigt stress vector.
!--------------------------------------------------------------------------------------------------
pure function math_33toVoigt6_stress(sigma) result(sigma_tilde)
real(pReal), dimension(6) :: sigma_tilde
real(pReal), dimension(3,3), intent(in) :: sigma
sigma_tilde = [sigma(1,1), sigma(2,2), sigma(3,3), &
sigma(3,2), sigma(3,1), sigma(1,2)]
end function math_33toVoigt6_stress
!--------------------------------------------------------------------------------------------------
!> @brief Convert 3x3 tensor into 6 Voigt strain vector.
!--------------------------------------------------------------------------------------------------
pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde)
real(pReal), dimension(6) :: epsilon_tilde
real(pReal), dimension(3,3), intent(in) :: epsilon
epsilon_tilde = [ epsilon(1,1), epsilon(2,2), epsilon(3,3), &
2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)]
end function math_33toVoigt6_strain
!--------------------------------------------------------------------------------------------------
!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix.
!--------------------------------------------------------------------------------------------------
pure function math_Voigt66to3333(m66)
@ -875,7 +938,7 @@ end function math_Voigt66to3333
!--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix
!> @brief Convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix.
!--------------------------------------------------------------------------------------------------
pure function math_3333toVoigt66(m3333)
@ -924,7 +987,7 @@ real(pReal) function math_sampleGaussVar(mu, sigma, width)
do
call random_number(rnd)
scatter = width_ * (2.0_pReal * rnd(1) - 1.0_pReal)
if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn
if (rnd(2) <= exp(-0.5_pReal * scatter**2)) exit ! test if scattered value is drawn
enddo
math_sampleGaussVar = scatter * sigma
@ -1086,15 +1149,15 @@ function math_eigvalsh33(m)
I = math_invariantsSym33(m) ! invariants are coefficients in characteristic polynomial apart for the sign of c0 and c2 in http://arxiv.org/abs/physics/0610206
P = I(2)-I(1)**2.0_pReal/3.0_pReal ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
P = I(2)-I(1)**2/3.0_pReal ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
Q = product(I(1:2))/3.0_pReal &
- 2.0_pReal/27.0_pReal*I(1)**3.0_pReal &
- 2.0_pReal/27.0_pReal*I(1)**3 &
- I(3) ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
if (all(abs([P,Q]) < TOL)) then
math_eigvalsh33 = math_eigvalsh(m)
else
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
rho=sqrt(-3.0_pReal*P**3)/9.0_pReal
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), &

View File

@ -5,6 +5,7 @@
!--------------------------------------------------------------------------------------------------
module phase
use prec
use constants
use math
use rotations
use IO

View File

@ -388,8 +388,7 @@ module function phase_K_phi(co,ce) result(K)
real(pReal), dimension(3,3) :: K
real(pReal), parameter :: l = 1.0_pReal
K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%D) &
* l**2.0_pReal
K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%D) * l**2
end function phase_K_phi

View File

@ -129,7 +129,7 @@ module subroutine anisobrittle_dotState(S, ph,en)
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal
traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2
damageState(ph)%dotState(1,en) = damageState(ph)%dotState(1,en) &
+ prm%dot_o / prm%s_crit(i) &
@ -190,7 +190,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
dLd_dTstar = 0.0_pReal
associate(prm => param(ph))
do i = 1,prm%sum_N_cl
traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal
traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2
traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
if (abs(traction_d) > traction_crit + tol_math_check) then

View File

@ -58,7 +58,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
associate(prm => param(kinematics_thermal_expansion_instance(p)))
kinematic_type => kinematics%get(k)
prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=0.0_pReal)
prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=T_ROOM)
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)
@ -98,14 +98,14 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
associate(prm => param(kinematics_thermal_expansion_instance(ph)))
Li = dot_T * ( &
prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
prm%A(1:3,1:3,1) & ! constant coefficient
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! 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 / 1. &
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
+ prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 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 &
)
end associate
dLi_dTstar = 0.0_pReal

View File

@ -1,13 +1,15 @@
submodule(phase:mechanical) elastic
type :: tParameters
real(pReal) :: &
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
end type tParameters
type(tParameters), allocatable, dimension(:) :: param
@ -28,7 +30,7 @@ module subroutine elastic_init(phases)
phase, &
mech, &
elastic
logical :: thermal_active
print'(/,1x,a)', '<<<+- phase:mechanical:elastic init -+>>>'
print'(/,1x,a)', '<<<+- phase:mechanical:elastic:Hooke init -+>>>'
@ -45,15 +47,35 @@ module subroutine elastic_init(phases)
associate(prm => param(ph))
prm%C_11 = elastic%get_asFloat('C_11')
prm%C_12 = elastic%get_asFloat('C_12')
prm%C_44 = elastic%get_asFloat('C_44')
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)
if (any(phase_lattice(ph) == ['hP','tI'])) then
prm%C_13 = elastic%get_asFloat('C_13')
prm%C_33 = elastic%get_asFloat('C_33')
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)
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 = elastic%get_asFloat('C_66')
end associate
end do
@ -69,21 +91,44 @@ module function elastic_C66(ph,en) result(C66)
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(6,6) :: C66
real(pReal) :: T
associate(prm => param(ph))
C66 = 0.0_pReal
C66(1,1) = prm%C_11
C66(1,2) = prm%C_12
C66(4,4) = prm%C_44
T = thermal_T(ph,en)
C66(1,1) = prm%C_11(1) &
+ prm%C_11(2)*(T - prm%T_ref)**1 &
+ prm%C_11(3)*(T - prm%T_ref)**2
C66(1,2) = prm%C_12(1) &
+ prm%C_12(2)*(T - prm%T_ref)**1 &
+ prm%C_12(3)*(T - prm%T_ref)**2
C66(4,4) = prm%C_44(1) &
+ prm%C_44(2)*(T - prm%T_ref)**1 &
+ prm%C_44(3)*(T - prm%T_ref)**2
if (any(phase_lattice(ph) == ['hP','tI'])) then
C66(1,3) = prm%C_13
C66(3,3) = prm%C_33
C66(1,3) = prm%C_13(1) &
+ prm%C_13(2)*(T - prm%T_ref)**1 &
+ prm%C_13(3)*(T - prm%T_ref)**2
C66(3,3) = prm%C_33(1) &
+ prm%C_33(2)*(T - prm%T_ref)**1 &
+ prm%C_33(3)*(T - prm%T_ref)**2
end if
if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66
if (phase_lattice(ph) == 'tI') then
C66(6,6) = prm%C_66(1) &
+ prm%C_66(2)*(T - prm%T_ref)**1 &
+ prm%C_66(3)*(T - prm%T_ref)**2
end if
C66 = lattice_symmetrize_C66(C66,phase_lattice(ph))
@ -148,15 +193,17 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E
real(pReal), dimension(6,6) :: C66
real(pReal), dimension(3,3,3,3) :: C
integer :: &
i, j
C = math_Voigt66to3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en))
C66 = phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)
C = math_Voigt66to3333(C66)
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice 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
do i =1,3; do j=1,3
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko

View File

@ -7,9 +7,6 @@
!--------------------------------------------------------------------------------------------------
submodule(phase:plastic) dislotungsten
real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
type :: tParameters
real(pReal) :: &
D = 1.0_pReal, & !< grain size
@ -169,7 +166,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%D = pl%get_asFloat('D')
prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl')
prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.)
@ -344,7 +341,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, &
0.0_pReal, &
prm%dipoleformation)
v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) &
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)) &
* (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
@ -475,7 +472,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
if (present(tau_pos_out)) tau_pos_out = tau_pos
if (present(tau_neg_out)) tau_neg_out = tau_neg
associate(BoltzmannRatio => prm%Q_s/(kB*T), &
associate(BoltzmannRatio => prm%Q_s/(K_B*T), &
b_rho_half => stt%rho_mob(:,en) * prm%b_sl * 0.5_pReal, &
effectiveLength => dst%Lambda_sl(:,en) - prm%w)
@ -501,7 +498,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
* StressRatio_pminus1 / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_pos
dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal
dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2
ddot_gamma_dtau_pos = b_rho_half * dvel
else where significantPositiveTau2
@ -531,7 +528,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
* StressRatio_pminus1 / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_neg
dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal
dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2
ddot_gamma_dtau_neg = b_rho_half * dvel
else where significantNegativeTau2

View File

@ -9,9 +9,6 @@
!--------------------------------------------------------------------------------------------------
submodule(phase:plastic) dislotwin
real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
type :: tParameters
real(pReal) :: &
Q_cl = 1.0_pReal, & !< activation energy for dislocation climb
@ -31,7 +28,7 @@ submodule(phase:plastic) dislotwin
delta_G = 1.0_pReal, & !< Free energy difference between austensite and martensite
i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation
h = 1.0_pReal, & !< Stack height of hex nucleus
T_ref = 0.0_pReal, &
T_ref = T_ROOM, &
a_cI = 1.0_pReal, &
a_cF = 1.0_pReal
real(pReal), dimension(2) :: &
@ -597,7 +594,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
shearBandingContribution: if (dNeq0(prm%v_sb)) then
E_kB_T = prm%E_sb/(kB*T)
E_kB_T = prm%E_sb/(K_B*T)
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
do i = 1,6
@ -694,8 +691,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
* (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (mu*prm%b_sl(i)), &
1.0_pReal, &
prm%ExtendedDislocations)
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
v_cl = 2.0_pReal*prm%omega*b_d**2*exp(-prm%Q_cl/(K_B*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3/(K_B*T)) - 1.0_pReal)
dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) &
/ (d_hat-prm%d_caron(i))
@ -787,14 +784,14 @@ module subroutine dislotwin_dependentState(T,ph,en)
+ 3.0_pReal*prm%b_tr*mu/(prm%L_tr*prm%b_tr) &
+ prm%h*prm%delta_G/(3.0_pReal*prm%b_tr)
dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2
dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2
x0 = mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu) ! ToDo: In the paper, this is the Burgers vector for slip
x0 = mu*prm%b_tw**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tw(:,en) = mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu) ! ToDo: In the paper, this is the Burgers vector for slip
x0 = mu*prm%b_tr**2/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tr(:,en) = mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
end associate
@ -907,7 +904,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
significantStress: where(tau_eff > tol_math_check)
stressRatio = tau_eff/prm%tau_0
StressRatio_p = stressRatio** prm%p
Q_kB_T = prm%Q_sl/(kB*T)
Q_kB_T = prm%Q_sl/(K_B*T)
v_wait_inverse = exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) &
/ prm%v_0
v_run_inverse = prm%B/(tau_eff*prm%b_sl)
@ -920,7 +917,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
/ prm%tau_0
dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff
dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) &
/ (v_wait_inverse+v_run_inverse)**2.0_pReal
/ (v_wait_inverse+v_run_inverse)**2
ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl
else where significantStress
dot_gamma_sl = 0.0_pReal
@ -980,7 +977,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
(prm%L_tw*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i))))
(1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tw(i,en)-tau(i))))
else
Ndot0=0.0_pReal
end if
@ -1049,7 +1046,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
(prm%L_tr*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i))))
(1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tr(i,en)-tau(i))))
else
Ndot0=0.0_pReal
end if

View File

@ -19,9 +19,6 @@ submodule(phase:plastic) nonlocal
type(tGeometry), dimension(:), allocatable :: geom
real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
! storage order of dislocation types
integer, dimension(*), parameter :: &
sgl = [1,2,3,4,5,6,7,8] !< signed (single)
@ -623,7 +620,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
* spread(( 1.0_pReal - prm%f_F &
+ prm%f_F &
* log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,en),prm%rho_significant))) &
/ log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl)
/ log(0.35_pReal * prm%b_sl * 1e6_pReal))**2,2,prm%sum_N_sl)
else
myInteractionMatrix = prm%h_sl_sl
end if
@ -1094,9 +1091,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
! thermally activated annihilation of edge dipoles by climb
rhoDotThermalAnnihilation = 0.0_pReal
D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) ! eq. 3.53
D_SD = prm%D_0 * exp(-prm%Q_cl / (K_B * Temperature)) ! eq. 3.53
v_climb = D_SD * mu * prm%V_at &
/ (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature) ! eq. 3.54
/ (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * K_B * Temperature) ! eq. 3.54
forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) &
rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
@ -1307,10 +1304,10 @@ function rhoDotFlux(timestep,ph,en,ip,el)
* math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (compatibility(c,:,s,n,ip,el) > 0.0_pReal) &
rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) &
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to equally signed mobile dislocation type
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to equally signed mobile dislocation type
where (compatibility(c,:,s,n,ip,el) < 0.0_pReal) &
rhoDotFlux(:,topp) = rhoDotFlux(:,topp) &
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to opposite signed mobile dislocation type
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to opposite signed mobile dislocation type
end if
end do
@ -1339,7 +1336,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
c = (t + 1) / 2
if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density
transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
transmissivity = sum(compatibility(c,:,s,n,ip,el)**2) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal
end if
@ -1666,14 +1663,14 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T,
!* Peierls contribution
tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s))
lambda_P = prm%b_sl(s)
activationVolume_P = prm%w *prm%b_sl(s)**3.0_pReal
activationVolume_P = prm%w *prm%b_sl(s)**3
criticalStress_P = prm%peierlsStress(s,c)
activationEnergy_P = criticalStress_P * activationVolume_P
tauRel_P = min(1.0_pReal, tauEff / criticalStress_P)
tPeierls = 1.0_pReal / prm%nu_a &
* exp(activationEnergy_P / (kB * T) &
* exp(activationEnergy_P / (K_B * T) &
* (1.0_pReal - tauRel_P**prm%p)**prm%q)
dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (kB * T) &
dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) &
* (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), &
0.0_pReal, &
tauEff < criticalStress_P)
@ -1681,12 +1678,12 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T,
! Contribution from solid solution strengthening
tauEff = abs(tau(s)) - tauThreshold(s)
lambda_S = prm%b_sl(s) / sqrt(prm%c_sol)
activationVolume_S = prm%f_sol * prm%b_sl(s)**3.0_pReal / sqrt(prm%c_sol)
activationVolume_S = prm%f_sol * prm%b_sl(s)**3 / sqrt(prm%c_sol)
criticalStress_S = prm%Q_sol / activationVolume_S
tauRel_S = min(1.0_pReal, tauEff / criticalStress_S)
tSolidSolution = 1.0_pReal / prm%nu_a &
* exp(prm%Q_sol / (kB * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * T) &
* exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) &
* (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), &
0.0_pReal, &
tauEff < criticalStress_S)
@ -1697,8 +1694,8 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T,
v(s) = sign(1.0_pReal,tau(s)) &
/ (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff))
dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2.0_pReal))
dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / lambda_P
dv_dtau(s) = v(s)**2 * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2))
dv_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P
end if
end do

View File

@ -75,6 +75,7 @@ module rotations
procedure, public :: rotVector
procedure, public :: rotTensor2
procedure, public :: rotTensor4
procedure, public :: rotStiffness
procedure, public :: misorientation
procedure, public :: standardize
end type rotation
@ -339,8 +340,7 @@ end function rotTensor2
!---------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief rotate a rank-4 tensor passively (default) or actively
!> @brief Rotate a rank-4 tensor passively (default) or actively.
!> @details: rotation is based on rotation matrix
!! ToDo: Need to check active/passive !!!
!---------------------------------------------------------------------------------------------------
@ -354,6 +354,7 @@ pure function rotTensor4(self,T,active) result(tRot)
real(pReal), dimension(3,3) :: R
integer :: i,j,k,l,m,n,o,p
if (present(active)) then
R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
else
@ -371,7 +372,47 @@ end function rotTensor4
!---------------------------------------------------------------------------------------------------
!> @brief misorientation
!> @brief Rotate a rank-4 tensor in Voigt 6x6 notation passively (default) or actively.
!> @details: https://scicomp.stackexchange.com/questions/35600
!! ToDo: Need to check active/passive !!!
!---------------------------------------------------------------------------------------------------
pure function rotStiffness(self,C,active) result(cRot)
real(pReal), dimension(6,6) :: cRot
class(rotation), intent(in) :: self
real(pReal), intent(in), dimension(6,6) :: C
logical, intent(in), optional :: active
real(pReal), dimension(3,3) :: R
real(pReal), dimension(6,6) :: M
if (present(active)) then
R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
else
R = self%asMatrix()
endif
M = reshape([R(1,1)**2.0_pReal, R(2,1)**2.0_pReal, R(3,1)**2.0_pReal, &
R(2,1)*R(3,1), R(1,1)*R(3,1), R(1,1)*R(2,1), &
R(1,2)**2.0_pReal, R(2,2)**2.0_pReal, R(3,2)**2.0_pReal, &
R(2,2)*R(3,2), R(1,2)*R(3,2), R(1,2)*R(2,2), &
R(1,3)**2.0_pReal, R(2,3)**2.0_pReal, R(3,3)**2.0_pReal, &
R(2,3)*R(3,3), R(1,3)*R(3,3), R(1,3)*R(2,3), &
2.0_pReal*R(1,2)*R(1,3), 2.0_pReal*R(2,2)*R(2,3), 2.0_pReal*R(3,2)*R(3,3), &
R(2,2)*R(3,3)+R(2,3)*R(3,2), R(1,2)*R(3,3)+R(1,3)*R(3,2), R(1,2)*R(2,3)+R(1,3)*R(2,2), &
2.0_pReal*R(1,3)*R(1,1), 2.0_pReal*R(2,3)*R(2,1), 2.0_pReal*R(3,3)*R(3,1), &
R(2,3)*R(3,1)+R(2,1)*R(3,3), R(1,3)*R(3,1)+R(1,1)*R(3,3), R(1,3)*R(2,1)+R(1,1)*R(2,3), &
2.0_pReal*R(1,1)*R(1,2), 2.0_pReal*R(2,1)*R(2,2), 2.0_pReal*R(3,1)*R(3,2), &
R(2,1)*R(3,2)+R(2,2)*R(3,1), R(1,1)*R(3,2)+R(1,2)*R(3,1), R(1,1)*R(2,2)+R(1,2)*R(2,1)],[6,6])
cRot = matmul(M,matmul(C,transpose(M)))
end function rotStiffness
!---------------------------------------------------------------------------------------------------
!> @brief Misorientation.
!---------------------------------------------------------------------------------------------------
pure elemental function misorientation(self,other)
@ -386,7 +427,7 @@ end function misorientation
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to rotation matrix
!> @brief Convert unit quaternion to rotation matrix.
!---------------------------------------------------------------------------------------------------
pure function qu2om(qu) result(om)
@ -395,8 +436,8 @@ pure function qu2om(qu) result(om)
real(pReal) :: qq
qq = qu(1)**2-sum(qu(2:4)**2)
qq = qu(1)**2-sum(qu(2:4)**2)
om(1,1) = qq+2.0_pReal*qu(2)**2
om(2,2) = qq+2.0_pReal*qu(3)**2
@ -416,7 +457,7 @@ end function qu2om
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to Euler angles
!> @brief Convert unit quaternion to Euler angles.
!---------------------------------------------------------------------------------------------------
pure function qu2eu(qu) result(eu)
@ -425,6 +466,7 @@ pure function qu2eu(qu) result(eu)
real(pReal) :: q12, q03, chi
q03 = qu(1)**2+qu(4)**2
q12 = qu(2)**2+qu(3)**2
chi = sqrt(q03*q12)
@ -578,7 +620,7 @@ pure function om2eu(om) result(eu)
real(pReal) :: zeta
if (dNeq(abs(om(3,3)),1.0_pReal,1.e-8_pReal)) then
zeta = 1.0_pReal/sqrt(math_clip(1.0_pReal-om(3,3)**2.0_pReal,1e-64_pReal,1.0_pReal))
zeta = 1.0_pReal/sqrt(math_clip(1.0_pReal-om(3,3)**2,1e-64_pReal,1.0_pReal))
eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), &
acos(math_clip(om(3,3),-1.0_pReal,1.0_pReal)), &
atan2(om(1,3)*zeta, om(2,3)*zeta)]
@ -1099,7 +1141,7 @@ pure function ho2ax(ho) result(ax)
+0.000003953714684212874_pReal, -0.00000036555001439719544_pReal ]
! normalize h and store the magnitude
hmag_squared = sum(ho**2.0_pReal)
hmag_squared = sum(ho**2)
if (dEq0(hmag_squared)) then
ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ]
else
@ -1379,6 +1421,7 @@ subroutine selfTest()
real(pReal), dimension(3) :: x, eu, ho, v3
real(pReal), dimension(3,3) :: om, t33
real(pReal), dimension(3,3,3,3) :: t3333
real(pReal), dimension(6,6) :: C
real :: A,B
integer :: i
@ -1412,6 +1455,7 @@ subroutine selfTest()
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu/qu2om'
if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu/qu2eu'
if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu/qu2ax'
@ -1447,17 +1491,22 @@ subroutine selfTest()
call R%fromMatrix(om)
call random_number(v3)
if(all(dNeq(R%rotVector(R%rotVector(v3),active=.true.),v3,1.0e-12_pReal))) &
if (any(dNeq(R%rotVector(R%rotVector(v3),active=.true.),v3,1.0e-12_pReal))) &
error stop 'rotVector'
call random_number(t33)
if(all(dNeq(R%rotTensor2(R%rotTensor2(t33),active=.true.),t33,1.0e-12_pReal))) &
if (any(dNeq(R%rotTensor2(R%rotTensor2(t33),active=.true.),t33,1.0e-12_pReal))) &
error stop 'rotTensor2'
call random_number(t3333)
if(all(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) &
if (any(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) &
error stop 'rotTensor4'
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))) &
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
end do