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 internal-restructure

This commit is contained in:
Sharan Roongta 2021-03-11 19:36:56 +01:00
parent 73d50c9bed 8daa8a965c
commit 65cd945aea
40 changed files with 533 additions and 534 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 0289c1bbfec1a1aef77a8cbaeed134035549e738 Subproject commit 13dfa0ee9d702782f0b7999f3f7fb2384f58d768

2
VERSION
View File

@ -1 +1 @@
v3.0.0-alpha2-534-g1d6234181 v3.0.0-alpha2-602-ge2d4ab427

4
examples/ConfigFiles/Microstructure_DP_Steel.config
View File

@ -1,4 +0,0 @@
[DP_Steel]
crystallite 1
(constituent) phase 1 texture 1 fraction 0.82
(constituent) phase 2 texture 2 fraction 0.18

64
examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config
View File

@ -1,64 +0,0 @@
[TWIP_Steel_FeMnC]
elasticity hooke
plasticity dislotwin
(output) rho_mob
(output) rho_dip
(output) gamma_sl
(output) lambda_sl
(output) tau_pass
(output) f_tw
(output) lambda_tw
(output) tau_hat_tw
(output) f_tr
### Material parameters ###
lattice_structure fcc
C11 175.0e9 # From Music et al. Applied Physics Letters 91, 191904 (2007)
C12 115.0e9
C44 135.0e9
grainsize 2.0e-5 # Average grain size [m]
SolidSolutionStrength 1.5e8 # Strength due to elements in solid solution
### Dislocation glide parameters ###
Nslip 12
slipburgers 2.56e-10 # Burgers vector of slip system [m]
rhoedgedip0 1.0 # Initial dislocation density [m/m**3]
rhoedge0 1.0e12 # Initial dislocation density [m/m**3]
v0 1.0e-4 # Initial glide velocity [m/s]
Qedge 3.7e-19 # Activation energy for dislocation glide [J]
p_slip 1.0 # p-exponent in glide velocity
q_slip 1.0 # q-exponent in glide velocity
# hardening of glide
CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path
D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s]
Qsd 4.5e-19 # Activation energy for climb [J]
Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b^3]
Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b]
interactionSlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Interaction coefficients (Kubin et al. 2008)
### Shearband parameters ###
shearbandresistance 180e6
shearbandvelocity 0e-4 # set to zero to turn shear banding of
QedgePerSbSystem 3.7e-19 # Activation energy for shear banding [J]
p_shearband 1.0 # p-exponent in glide velocity
q_shearband 1.0 # q-exponent in glide velocity
### Twinning parameters ###
Ntwin 12
twinburgers 1.47e-10 # Burgers vector of twin system [m]
twinsize 5.0e-8 # Twin stack mean thickness [m]
L0_twin 442.0 # Length of twin nuclei in Burgers vectors
maxtwinfraction 1.0 # Maximum admissible twin volume fraction
xc_twin 1.0e-9 # critical distance for formation of twin nucleus
VcrossSlip 1.67e-29 # cross slip volume
r_twin 10.0 # r-exponent in twin formation probability
Cmfptwin 1.0 # Adj. parameter controlling twin mean free path
Cthresholdtwin 1.0 # Adj. parameter controlling twin threshold stress
interactionSlipTwin 0.0 1.0 1.0 # Dislocation-Twin interaction coefficients
interactionTwinTwin 0.0 1.0 # Twin-Twin interaction coefficients
SFE_0K -0.0396 # stacking fault energy at zero K; TWIP steel: -0.0526; Cu: -0.0396
dSFE_dT 0.0002 # temperature dependance of stacking fault energy

41
examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml Normal file
View File

@ -0,0 +1,41 @@
TWIP_Steel_FeMnC:
lattice: cF
mechanics:
elasticity: {type: hooke, C_11: 175.0e9, C_12: 115.0e9, C_44: 135.0e9}
plasticity:
type: dislotwin
output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, tau_hat_tw, f_tr]
D: 2.0e-5
N_sl: [12]
b_sl: [2.56e-10]
rho_mob_0: [1.0e12]
rho_dip_0: [1.0]
v_0: [1.0e4]
Q_s: [3.7e-19]
p_sl: [1.0]
q_sl: [1.0]
tau_0: [1.5e8]
i_sl: [10.0] # Adj. parameter controlling dislocation mean free path
D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s
D_a: 1.0 # minimum dipole distance / b
Q_cl: 4.5e-19 # Activation energy for climb / J
h_sl_sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.122] # Interaction coefficients (Kubin et al. 2008)
# shear band parameters
xi_sb: 180.0e6
Q_sb: 3.7e-19
p_sb: 1.0
q_sb: 1.0
v_sb: 0.0 # set to 0, to turn it off
# twinning parameters
N_tw: [12]
b_tw: [1.47e-10] # Burgers vector length of twin system / b
t_tw: [5.0e-8] # Twin stack mean thickness / m
L_tw: 442.0 # Length of twin nuclei / b
x_c_tw: 1.0e-9 # critical distance for formation of twin nucleus / m
V_cs: 1.67e-29 # cross slip volume / m^3
p_tw: [10.0] # r-exponent in twin formation probability
i_tw: 1.0 # Adj. parameter controlling twin mean free path
h_sl_tw: [0.0, 1.0, 1.0] # dislocation-twin interaction coefficients
h_tw_tw: [0.0, 1.0] # twin-twin interaction coefficients
Gamma_sf_0K: -0.0396 # stacking fault energy / J/m^2 at zero K; TWIP steel: -0.0526; Cu: -0.0396
dGamma_sf_dT: 0.0002 # temperature dependence / J/(m^2 K) of stacking fault energy

36
examples/ConfigFiles/Phase_Dislotwin_Tungsten.config
View File

@ -1,36 +0,0 @@
[Tungsten]
elasticity hooke
plasticity dislotwin
### Material parameters ###
lattice_structure bcc
C11 523.0e9 # From Marinica et al. Journal of Physics: Condensed Matter(2013)
C12 202.0e9
C44 161.0e9
grainsize 2.0e-5 # Average grain size [m]
SolidSolutionStrength 1.5e8 # Strength due to elements in solid solution
### Dislocation glide parameters ###
#per family
Nslip 12
slipburgers 2.72e-10 # Burgers vector of slip system [m]
rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3]
rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3]
v0 1.0e-4 # Initial glide velocity [m/s]
Qedge 2.725e-19 # Activation energy for dislocation glide [J]
p_slip 0.78 # p-exponent in glide velocity
q_slip 1.58 # q-exponent in glide velocity
tau_peierls 2.03e9 # peierls stress (for bcc)
dipoleformationfactor 0 # to have hardening due to dipole formation off
#hardening
CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path
D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s]
Qsd 4.5e-19 # Activation energy for climb [J]
Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b]
Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b]
interaction_slipslip 1 1 1.4 1.4 1.4 1.4

21
examples/ConfigFiles/Phase_Dislotwin_Tungsten.yaml Normal file
View File

@ -0,0 +1,21 @@
Tungsten:
lattice: cI
mechanics:
elasticity: {type: hooke, C_11: 523.0e9, C_12: 202.0e9, C_44: 161.0e9} # Marinica et al. Journal of Physics: Condensed Matter(2013)
plasticity:
type: dislotwin
D: 2.0e-5 # Average grain size / m
N_sl: [12]
b_sl: [2.72e-10] # Burgers vector length of slip families / m
rho_mob_0: [1.0e12]
rho_dip_0: [1.0]
v_0: [1.0e4] # Initial glide velocity / m/s
Q_s: [2.725e-19] # Activation energy for dislocation glide / J
p_sl: [0.78] # p-exponent in glide velocity
q_sl: [1.58] # q-exponent in glide velocity
tau_0: [1.5e8] # solid solution strength / Pa
i_sl: [10.0] # Adj. parameter controlling dislocation mean free path
D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s
D_a: 1.0 # minimum dipole distance / b
Q_cl: 4.5e-19 # Activation energy for climb / J
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]

3
examples/ConfigFiles/Phase_Hydrogen.config
View File

@ -1,3 +0,0 @@
hydrogenflux_diffusion11 1.0
hydrogenflux_mobility11 1.0
hydrogenVolume 1e-28

9
examples/ConfigFiles/Phase_Phenopowerlaw_cpTi.yaml
View File

@ -1,13 +1,14 @@
# M. Levy, Handbook of Elastic Properties of Solids, Liquids, and Gases (2001) # M. Levy, Handbook of Elastic Properties of Solids, Liquids, and Gases (2001)
# C. Zambaldi, "Orientation informed nanoindentation of a-titanium: Indentation pileup in hexagonal metals deforming by prismatic slip", J. Mater. Res., Vol. 27, No. 1, Jan 14, 2012 # C. Zambaldi, "Orientation informed nanoindentation of a-titanium: Indentation pileup in hexagonal metals deforming by prismatic slip", J. Mater. Res., Vol. 27, No. 1, Jan 14, 2012
Ti-alpha: # Better use values from L. Wang, Z. Zheng, H. Phukan, P. Kenesei, J.-S. Park, J. Lind, R.M. Suter, T.R. Bieler, Direct measurement of critical resolved shear stress of prismatic and basal slip in polycrystalline Ti using high energy X-ray diffraction microscopy, Acta Mater 2017
cpTi:
lattice: hP lattice: hP
c/a: 1.587 c/a: 1.587
mechanics: mechanics:
output: [F, P, F_e, F_p, L_p, O] output: [F, P, F_e, F_p, L_p, O]
elasticity: {C_11: 160.0e9, C_12: 90.0e9, C_13: 66.0e9, C_33: 181.7e9, C_44: 46.5e9, type: hooke} elasticity: {C_11: 160.0e9, C_12: 90.0e9, C_13: 66.0e9, C_33: 181.7e9, C_44: 46.5e9, type: hooke}
plasticity: plasticity:
N_sl: [3, 3, 0, 0, 12] N_sl: [3, 3, 0, 6, 12]
a_sl: 2.0 a_sl: 2.0
dot_gamma_0_sl: 0.001 dot_gamma_0_sl: 0.001
h_0_sl_sl: 200e6 h_0_sl_sl: 200e6
@ -15,5 +16,5 @@ Ti-alpha:
n_sl: 20 n_sl: 20
output: [gamma_sl] output: [gamma_sl]
type: phenopowerlaw type: phenopowerlaw
xi_0_sl: [349e6, 150e6, 0, 0, 1107e6] xi_0_sl: [0.15e9, 0.09e9, 0, 0.20e9, 0.25e9]
xi_inf_sl: [568e6, 1502e6, 0, 0, 3420e6] xi_inf_sl: [0.24e9, 0.5e9, 0, 0.6e9, 0.8e9]

61
python/damask/_config.py
View File

@ -1,5 +1,6 @@
import copy import copy
from io import StringIO from io import StringIO
from collections.abc import Iterable
import abc import abc
import numpy as np import numpy as np
@ -44,6 +45,42 @@ class Config(dict):
copy = __copy__ copy = __copy__
def __or__(self,other):
"""
Update configuration with contents of other.
Parameters
----------
other : damask.Config or dict
Key-value pairs that update self.
"""
duplicate = self.copy()
duplicate.update(other)
return duplicate
def __ior__(self,other):
"""Update configuration with contents of other."""
return self.__or__(other)
def delete(self,keys):
"""
Remove configuration keys.
Parameters
----------
keys : iterable or scalar
Label of the key(s) to remove.
"""
duplicate = self.copy()
for k in keys if isinstance(keys, Iterable) and not isinstance(keys, str) else [keys]:
del duplicate[k]
return duplicate
@classmethod @classmethod
def load(cls,fname): def load(cls,fname):
""" """
@ -99,30 +136,6 @@ class Config(dict):
fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs)) fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
def add(self,d):
"""
Add dictionary.
d : dict
Dictionary to append.
"""
duplicate = self.copy()
duplicate.update(d)
return duplicate
def delete(self,key):
"""
Delete item.
key : str or scalar
Label of the key to remove.
"""
duplicate = self.copy()
del duplicate[key]
return duplicate
@property @property
@abc.abstractmethod @abc.abstractmethod
def is_complete(self): def is_complete(self):

8
python/damask/_result.py
View File

@ -371,7 +371,7 @@ class Result:
with h5py.File(self.fname,'r') as f: with h5py.File(self.fname,'r') as f:
for dataset in sets: for dataset in sets:
for group in self.groups_with_datasets(dataset): for group in self.groups_with_datasets(dataset):
path = os.path.join(group,dataset) path = '/'.join([group,dataset])
inc,prop,name,cat,item = (path.split('/') + ['']*5)[:5] inc,prop,name,cat,item = (path.split('/') + ['']*5)[:5]
key = '/'.join([prop,name+tag]) key = '/'.join([prop,name+tag])
if key not in inGeom: if key not in inGeom:
@ -388,15 +388,15 @@ class Result:
np.nan, np.nan,
dtype=np.dtype(f[path])) dtype=np.dtype(f[path]))
data[inGeom[key]] = (f[path] if len(shape)>1 else np.expand_dims(f[path],1))[inData[key]] data[inGeom[key]] = (f[path] if len(shape)>1 else np.expand_dims(f[path],1))[inData[key]]
path = (os.path.join(*([prop,name]+([cat] if cat else [])+([item] if item else []))) if split else path)+tag path = ('/'.join([prop,name]+([cat] if cat else [])+([item] if item else [])) if split else path)+tag
if split: if split:
try: try:
tbl[inc].add(path,data) tbl[inc] = tbl[inc].add(path,data)
except KeyError: except KeyError:
tbl[inc] = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]}) tbl[inc] = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]})
else: else:
try: try:
tbl.add(path,data) tbl = tbl.add(path,data)
except AttributeError: except AttributeError:
tbl = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]}) tbl = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]})

14
python/damask/_rotation.py
View File

@ -475,15 +475,17 @@ class Rotation:
Parameters Parameters
---------- ----------
vector : bool, optional compact : bool, optional
Return as actual Rodrigues-Frank vector, i.e. axis Return as actual Rodrigues-Frank vector,
and angle argument are not separated. i.e. axis and angle argument are not separated.
Returns Returns
------- -------
rho : numpy.ndarray of shape (...,4) unless vector == True: rho : numpy.ndarray of shape (...,4) containing
numpy.ndarray of shape (...,3) [n_1, n_2, n_3, tan(ω/2)], ǀnǀ = 1 and ω ∈ [0,π]
Rodrigues-Frank vector: [n_1, n_2, n_3, tan(ω/2)], ǀnǀ = 1 and ω ∈ [0,π]. unless compact == True:
numpy.ndarray of shape (...,3) containing
tan(ω/2) [n_1, n_2, n_3], ω ∈ [0,π].
""" """
ro = Rotation._qu2ro(self.quaternion) ro = Rotation._qu2ro(self.quaternion)

40
python/damask/_table.py
View File

@ -189,6 +189,11 @@ class Table:
label : str label : str
Column label. Column label.
Returns
-------
data : numpy.ndarray
Array of column data.
""" """
if re.match(r'[0-9]*?_',label): if re.match(r'[0-9]*?_',label):
idx,key = label.split('_',1) idx,key = label.split('_',1)
@ -212,6 +217,11 @@ class Table:
info : str, optional info : str, optional
Human-readable information about the new data. Human-readable information about the new data.
Returns
-------
table : Table
Updated table.
""" """
dup = self.copy() dup = self.copy()
dup._add_comment(label,data.shape[1:],info) dup._add_comment(label,data.shape[1:],info)
@ -238,6 +248,11 @@ class Table:
info : str, optional info : str, optional
Human-readable information about the modified data. Human-readable information about the modified data.
Returns
-------
table : Table
Updated table.
""" """
dup = self.copy() dup = self.copy()
dup._add_comment(label,data.shape[1:],info) dup._add_comment(label,data.shape[1:],info)
@ -261,6 +276,11 @@ class Table:
label : str label : str
Column label. Column label.
Returns
-------
table : Table
Updated table.
""" """
dup = self.copy() dup = self.copy()
dup.data.drop(columns=label,inplace=True) dup.data.drop(columns=label,inplace=True)
@ -279,6 +299,11 @@ class Table:
label_new : str or iterable of str label_new : str or iterable of str
New column label(s). New column label(s).
Returns
-------
table : Table
Updated table.
""" """
dup = self.copy() dup = self.copy()
columns = dict(zip([old] if isinstance(old,str) else old, columns = dict(zip([old] if isinstance(old,str) else old,
@ -300,6 +325,11 @@ class Table:
ascending : bool or list, optional ascending : bool or list, optional
Set sort order. Set sort order.
Returns
-------
table : Table
Updated table.
""" """
dup = self.copy() dup = self.copy()
dup._label_discrete() dup._label_discrete()
@ -320,6 +350,11 @@ class Table:
other : Table other : Table
Table to append. Table to append.
Returns
-------
table : Table
Concatenated table.
""" """
if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns): if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns):
raise KeyError('Labels or shapes or order do not match') raise KeyError('Labels or shapes or order do not match')
@ -340,6 +375,11 @@ class Table:
other : Table other : Table
Table to join. Table to join.
Returns
-------
table : Table
Joined table.
""" """
if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]: if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]:
raise KeyError('Dublicated keys or row count mismatch') raise KeyError('Dublicated keys or row count mismatch')

33
python/damask/util.py
View File

@ -133,6 +133,8 @@ def execute(cmd,
stdout = stdout.decode('utf-8').replace('\x08','') stdout = stdout.decode('utf-8').replace('\x08','')
stderr = stderr.decode('utf-8').replace('\x08','') stderr = stderr.decode('utf-8').replace('\x08','')
if process.returncode != 0: if process.returncode != 0:
print(stdout)
print(stderr)
raise RuntimeError(f"'{cmd}' failed with returncode {process.returncode}") raise RuntimeError(f"'{cmd}' failed with returncode {process.returncode}")
return stdout, stderr return stdout, stderr
@ -193,7 +195,7 @@ def scale_to_coprime(v):
return m return m
def project_stereographic(vector,normalize=False): def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
""" """
Apply stereographic projection to vector. Apply stereographic projection to vector.
@ -201,18 +203,37 @@ def project_stereographic(vector,normalize=False):
---------- ----------
vector : numpy.ndarray of shape (...,3) vector : numpy.ndarray of shape (...,3)
Vector coordinates to be projected. Vector coordinates to be projected.
direction : str
Projection direction 'x', 'y', or 'z'.
Defaults to 'z'.
normalize : bool normalize : bool
Ensure unit length for vector. Defaults to False. Ensure unit length of input vector. Defaults to True.
keepdims : bool
Maintain three-dimensional output coordinates.
Default 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.
Returns Returns
------- -------
coordinates : numpy.ndarray of shape (...,2) coordinates : numpy.ndarray of shape (...,2 | 3)
Projected coordinates. Projected coordinates.
Examples
--------
>>> project_stereographic(np.ones(3))
[0.3660254, 0.3660254]
>>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True)
[0, 0.5, 0.5]
>>> project_stereographic([0,1,1],direction='y',normalize=True,keepdims=False)
[0.41421356, 0]
""" """
v_ = vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector shift = 'zyx'.index(direction)
return np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])), v_ = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
np.zeros_like(v_[...,2:3])]) shift,axis=-1)
return np.roll(np.block([v_[...,:2]/(1+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,function_name=None): def execution_stamp(class_name,function_name=None):

11
python/tests/test_Config.py
View File

@ -23,8 +23,17 @@ class TestConfig:
assert Config.load(f) == config assert Config.load(f) == config
def test_add_remove(self): def test_add_remove(self):
dummy = {'hello':'world','foo':'bar'}
config = Config() config = Config()
assert config.add({'hello':'world'}).delete('hello') == config config |= dummy
assert config == Config() | dummy
config = config.delete(dummy)
assert config == Config()
assert (config | dummy ).delete( 'hello' ) == config | {'foo':'bar'}
assert (config | dummy ).delete([ 'hello', 'foo' ]) == config
assert (config | Config(dummy)).delete({ 'hello':1,'foo':2 }) == config
assert (config | Config(dummy)).delete(Config({'hello':1 })) == config | {'foo':'bar'}
def test_repr(self,tmp_path): def test_repr(self,tmp_path):
config = Config() config = Config()

19
python/tests/test_util.py
View File

@ -49,17 +49,18 @@ class TestUtil:
dist_sampled = np.histogram(centers[selected],bins)[0]/N_samples*np.sum(dist) dist_sampled = np.histogram(centers[selected],bins)[0]/N_samples*np.sum(dist)
assert np.sqrt(((dist - dist_sampled) ** 2).mean()) < .025 and selected.shape[0]==N_samples assert np.sqrt(((dist - dist_sampled) ** 2).mean()) < .025 and selected.shape[0]==N_samples
@pytest.mark.parametrize('point,normalize,answer', @pytest.mark.parametrize('point,direction,normalize,keepdims,answer',
[ [
([1,0,0],False,[1,0,0]), ([1,0,0],'z',False,True, [1,0,0]),
([1,0,0],True, [1,0,0]), ([1,0,0],'z',True, False,[1,0]),
([0,1,1],False,[0,0.5,0]), ([0,1,1],'z',False,True, [0,0.5,0]),
([0,1,1],True, [0,0.41421356,0]), ([0,1,1],'y',True, False,[0.41421356,0]),
([1,1,1],False,[0.5,0.5,0]), ([1,1,0],'x',False,False,[0.5,0]),
([1,1,1],True, [0.3660254, 0.3660254, 0]), ([1,1,1],'y',True, True, [0.3660254, 0,0.3660254]),
]) ])
def test_project_stereographic(self,point,normalize,answer): def test_project_stereographic(self,point,direction,normalize,keepdims,answer):
assert np.allclose(util.project_stereographic(np.array(point),normalize=normalize),answer) assert np.allclose(util.project_stereographic(np.array(point),direction=direction,
normalize=normalize,keepdims=keepdims),answer)
@pytest.mark.parametrize('fro,to,mode,answer', @pytest.mark.parametrize('fro,to,mode,answer',
[ [

6
src/CPFEM.f90
View File

@ -178,11 +178,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward
chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) !chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
!case (THERMAL_conduction_ID) chosenThermal1 ! case (THERMAL_conduction_ID) chosenThermal1
! temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & ! temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
! temperature_inp ! temperature_inp
end select chosenThermal1 !end select chosenThermal1
homogenization_F0(1:3,1:3,ma) = ffn homogenization_F0(1:3,1:3,ma) = ffn
homogenization_F(1:3,1:3,ma) = ffn1 homogenization_F(1:3,1:3,ma) = ffn1

100
src/YAML_types.f90
View File

@ -207,11 +207,11 @@ subroutine selfTest
select type(s1) select type(s1)
class is(tScalar) class is(tScalar)
s1 = '1' s1 = '1'
if(s1%asInt() /= 1) error stop 'tScalar_asInt' if (s1%asInt() /= 1) error stop 'tScalar_asInt'
if(dNeq(s1%asFloat(),1.0_pReal)) error stop 'tScalar_asFloat' if (dNeq(s1%asFloat(),1.0_pReal)) error stop 'tScalar_asFloat'
s1 = 'true' s1 = 'true'
if(.not. s1%asBool()) error stop 'tScalar_asBool' if (.not. s1%asBool()) error stop 'tScalar_asBool'
if(s1%asString() /= 'true') error stop 'tScalar_asString' if (s1%asString() /= 'true') error stop 'tScalar_asString'
end select end select
block block
@ -232,18 +232,18 @@ subroutine selfTest
call l1%append(s1) call l1%append(s1)
call l1%append(s2) call l1%append(s2)
n => l1 n => l1
if(any(l1%asInts() /= [2,3])) error stop 'tList_asInts' if (any(l1%asInts() /= [2,3])) error stop 'tList_asInts'
if(any(dNeq(l1%asFloats(),[2.0_pReal,3.0_pReal]))) error stop 'tList_asFloats' if (any(dNeq(l1%asFloats(),[2.0_pReal,3.0_pReal]))) error stop 'tList_asFloats'
if(n%get_asInt(1) /= 2) error stop 'byIndex_asInt' if (n%get_asInt(1) /= 2) error stop 'byIndex_asInt'
if(dNeq(n%get_asFloat(2),3.0_pReal)) error stop 'byIndex_asFloat' if (dNeq(n%get_asFloat(2),3.0_pReal)) error stop 'byIndex_asFloat'
endselect endselect
allocate(tList::l2) allocate(tList::l2)
select type(l2) select type(l2)
class is(tList) class is(tList)
call l2%append(l1) call l2%append(l1)
if(any(l2%get_asInts(1) /= [2,3])) error stop 'byIndex_asInts' if (any(l2%get_asInts(1) /= [2,3])) error stop 'byIndex_asInts'
if(any(dNeq(l2%get_asFloats(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_asFloats' if (any(dNeq(l2%get_asFloats(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_asFloats'
n => l2 n => l2
end select end select
deallocate(n) deallocate(n)
@ -265,10 +265,10 @@ subroutine selfTest
call l1%append(s2) call l1%append(s2)
n => l1 n => l1
if(any(l1%asBools() .neqv. [.true., .false.])) error stop 'tList_asBools' if (any(l1%asBools() .neqv. [.true., .false.])) error stop 'tList_asBools'
if(any(l1%asStrings() /= ['true ','False'])) error stop 'tList_asStrings' if (any(l1%asStrings() /= ['true ','False'])) error stop 'tList_asStrings'
if(n%get_asBool(2)) error stop 'byIndex_asBool' if (n%get_asBool(2)) error stop 'byIndex_asBool'
if(n%get_asString(1) /= 'true') error stop 'byIndex_asString' if (n%get_asString(1) /= 'true') error stop 'byIndex_asString'
end block end block
end subroutine selfTest end subroutine selfTest
@ -418,7 +418,7 @@ function tNode_get_byIndex(self,i) result(node)
integer :: j integer :: j
self_ => self%asList() self_ => self%asList()
if(i < 1 .or. i > self_%length) call IO_error(150,ext_msg='tNode_get_byIndex') if (i < 1 .or. i > self_%length) call IO_error(150,ext_msg='tNode_get_byIndex')
j = 1 j = 1
item => self_%first item => self_%first
@ -599,7 +599,7 @@ function tNode_getKey_byIndex(self,i) result(key)
dict => self%asDict() dict => self%asDict()
item => dict%first item => dict%first
do j = 1, dict%length do j = 1, dict%length
if(j == i) then if (j == i) then
key = item%key key = item%key
exit exit
else else
@ -624,18 +624,18 @@ function tNode_contains(self,k) result(exists)
type(tDict), pointer :: dict type(tDict), pointer :: dict
exists = .false. exists = .false.
if(self%isDict()) then if (self%isDict()) then
dict => self%asDict() dict => self%asDict()
do j=1, dict%length do j=1, dict%length
if(dict%getKey(j) == k) then if (dict%getKey(j) == k) then
exists = .true. exists = .true.
return return
endif endif
enddo enddo
elseif(self%isList()) then elseif (self%isList()) then
list => self%asList() list => self%asList()
do j =1, list%length do j=1, list%length
if(list%get_asString(j) == k) then if (list%get_asString(j) == k) then
exists = .true. exists = .true.
return return
endif endif
@ -663,7 +663,7 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
logical :: found logical :: found
found = present(defaultVal) found = present(defaultVal)
if(found) node => defaultVal if (found) node => defaultVal
self_ => self%asDict() self_ => self%asDict()
@ -681,7 +681,7 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
if (.not. found) then if (.not. found) then
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
else else
if(associated(item)) node => item%node if (associated(item)) node => item%node
endif endif
end function tNode_get_byKey end function tNode_get_byKey
@ -700,11 +700,11 @@ function tNode_get_byKey_asFloat(self,k,defaultVal) result(nodeAsFloat)
class(tNode), pointer :: node class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
if(self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsFloat = scalar%asFloat() nodeAsFloat = scalar%asFloat()
elseif(present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsFloat = defaultVal nodeAsFloat = defaultVal
else else
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
@ -726,11 +726,11 @@ function tNode_get_byKey_asInt(self,k,defaultVal) result(nodeAsInt)
class(tNode), pointer :: node class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
if(self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsInt = scalar%asInt() nodeAsInt = scalar%asInt()
elseif(present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsInt = defaultVal nodeAsInt = defaultVal
else else
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
@ -752,11 +752,11 @@ function tNode_get_byKey_asBool(self,k,defaultVal) result(nodeAsBool)
class(tNode), pointer :: node class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
if(self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsBool = scalar%asBool() nodeAsBool = scalar%asBool()
elseif(present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsBool = defaultVal nodeAsBool = defaultVal
else else
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
@ -778,11 +778,11 @@ function tNode_get_byKey_asString(self,k,defaultVal) result(nodeAsString)
class(tNode), pointer :: node class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
if(self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsString = scalar%asString() nodeAsString = scalar%asString()
elseif(present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsString = defaultVal nodeAsString = defaultVal
else else
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
@ -806,18 +806,18 @@ function tNode_get_byKey_asFloats(self,k,defaultVal,requiredSize) result(nodeAsF
class(tNode), pointer :: node class(tNode), pointer :: node
type(tList), pointer :: list type(tList), pointer :: list
if(self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
list => node%asList() list => node%asList()
nodeAsFloats = list%asFloats() nodeAsFloats = list%asFloats()
elseif(present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsFloats = defaultVal nodeAsFloats = defaultVal
else else
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
endif endif
if(present(requiredSize)) then if (present(requiredSize)) then
if(requiredSize /= size(nodeAsFloats)) call IO_error(146,ext_msg=k) if (requiredSize /= size(nodeAsFloats)) call IO_error(146,ext_msg=k)
endif endif
end function tNode_get_byKey_asFloats end function tNode_get_byKey_asFloats
@ -837,18 +837,18 @@ function tNode_get_byKey_asInts(self,k,defaultVal,requiredSize) result(nodeAsInt
class(tNode), pointer :: node class(tNode), pointer :: node
type(tList), pointer :: list type(tList), pointer :: list
if(self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
list => node%asList() list => node%asList()
nodeAsInts = list%asInts() nodeAsInts = list%asInts()
elseif(present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsInts = defaultVal nodeAsInts = defaultVal
else else
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
endif endif
if(present(requiredSize)) then if (present(requiredSize)) then
if(requiredSize /= size(nodeAsInts)) call IO_error(146,ext_msg=k) if (requiredSize /= size(nodeAsInts)) call IO_error(146,ext_msg=k)
endif endif
end function tNode_get_byKey_asInts end function tNode_get_byKey_asInts
@ -867,11 +867,11 @@ function tNode_get_byKey_asBools(self,k,defaultVal) result(nodeAsBools)
class(tNode), pointer :: node class(tNode), pointer :: node
type(tList), pointer :: list type(tList), pointer :: list
if(self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
list => node%asList() list => node%asList()
nodeAsBools = list%asBools() nodeAsBools = list%asBools()
elseif(present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsBools = defaultVal nodeAsBools = defaultVal
else else
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
@ -893,11 +893,11 @@ function tNode_get_byKey_asStrings(self,k,defaultVal) result(nodeAsStrings)
class(tNode), pointer :: node class(tNode), pointer :: node
type(tList), pointer :: list type(tList), pointer :: list
if(self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
list => node%asList() list => node%asList()
nodeAsStrings = list%asStrings() nodeAsStrings = list%asStrings()
elseif(present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsStrings = defaultVal nodeAsStrings = defaultVal
else else
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
@ -944,7 +944,7 @@ function tNode_get_byKey_asIndex(self,key) result(keyIndex)
item => dict%first item => dict%first
keyIndex = -1 keyIndex = -1
do i = 1, dict%length do i = 1, dict%length
if(key == item%key) then if (key == item%key) then
keyIndex = i keyIndex = i
exit exit
else else
@ -952,7 +952,7 @@ function tNode_get_byKey_asIndex(self,key) result(keyIndex)
endif endif
enddo enddo
if(keyIndex == -1) call IO_error(140,ext_msg=key) if (keyIndex == -1) call IO_error(140,ext_msg=key)
end function tNode_get_byKey_asIndex end function tNode_get_byKey_asIndex
@ -985,7 +985,7 @@ recursive function tList_asFormattedString(self,indent) result(str)
integer :: i, indent_ integer :: i, indent_
str = '' str = ''
if(present(indent)) then if (present(indent)) then
indent_ = indent indent_ = indent
else else
indent_ = 0 indent_ = 0
@ -993,7 +993,7 @@ recursive function tList_asFormattedString(self,indent) result(str)
item => self%first item => self%first
do i = 1, self%length do i = 1, self%length
if(i /= 1) str = str//repeat(' ',indent_) if (i /= 1) str = str//repeat(' ',indent_)
str = str//'- '//item%node%asFormattedString(indent_+2) str = str//'- '//item%node%asFormattedString(indent_+2)
item => item%next item => item%next
end do end do
@ -1014,7 +1014,7 @@ recursive function tDict_asFormattedString(self,indent) result(str)
integer :: i, indent_ integer :: i, indent_
str = '' str = ''
if(present(indent)) then if (present(indent)) then
indent_ = indent indent_ = indent
else else
indent_ = 0 indent_ = 0
@ -1022,7 +1022,7 @@ recursive function tDict_asFormattedString(self,indent) result(str)
item => self%first item => self%first
do i = 1, self%length do i = 1, self%length
if(i /= 1) str = str//repeat(' ',indent_) if (i /= 1) str = str//repeat(' ',indent_)
select type(node_1 =>item%node) select type(node_1 =>item%node)
class is(tScalar) class is(tScalar)
str = str//trim(item%key)//': '//item%node%asFormattedString(indent_+len_trim(item%key)+2) str = str//trim(item%key)//': '//item%node%asFormattedString(indent_+len_trim(item%key)+2)
@ -1270,7 +1270,7 @@ recursive subroutine tItem_finalize(self)
type(tItem),intent(inout) :: self type(tItem),intent(inout) :: self
deallocate(self%node) deallocate(self%node)
if(associated(self%next)) deallocate(self%next) if (associated(self%next)) deallocate(self%next)
end subroutine tItem_finalize end subroutine tItem_finalize

87
src/grid/DAMASK_grid.f90
View File

@ -37,9 +37,10 @@ program DAMASK_grid
f_out, & !< frequency of result writes f_out, & !< frequency of result writes
f_restart !< frequency of restart writes f_restart !< frequency of restart writes
logical :: estimate_rate !< follow trajectory of former loadcase logical :: estimate_rate !< follow trajectory of former loadcase
integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
end type tLoadCase end type tLoadCase
integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file ! variables related to information from load case and geom file
real(pReal), dimension(9) :: temp_valueVector !< temporarily from loadcase file when reading in tensors (initialize to 0.0) real(pReal), dimension(9) :: temp_valueVector !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
@ -53,6 +54,7 @@ program DAMASK_grid
integer, parameter :: & integer, parameter :: &
subStepFactor = 2 !< for each substep, divide the last time increment by 2.0 subStepFactor = 2 !< for each substep, divide the last time increment by 2.0
real(pReal) :: & real(pReal) :: &
T_0 = 300.0_pReal, &
time = 0.0_pReal, & !< elapsed time time = 0.0_pReal, & !< elapsed time
time0 = 0.0_pReal, & !< begin of interval time0 = 0.0_pReal, & !< begin of interval
timeinc = 1.0_pReal, & !< current time interval timeinc = 1.0_pReal, & !< current time interval
@ -78,8 +80,7 @@ program DAMASK_grid
maxCutBack, & !< max number of cut backs maxCutBack, & !< max number of cut backs
stagItMax !< max number of field level staggered iterations stagItMax !< max number of field level staggered iterations
character(len=pStringLen) :: & character(len=pStringLen) :: &
incInfo, & incInfo
loadcase_string
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
@ -98,10 +99,13 @@ program DAMASK_grid
quit quit
class (tNode), pointer :: & class (tNode), pointer :: &
num_grid, & num_grid, &
debug_grid, & ! pointer to grid debug options
config_load, & config_load, &
load_steps, & load_steps, &
load_step, & load_step, &
solver, &
initial_conditions, &
ic_thermal, &
thermal, &
step_bc, & step_bc, &
step_mech, & step_mech, &
step_discretization, & step_discretization, &
@ -112,17 +116,11 @@ program DAMASK_grid
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll call CPFEM_initAll
print'(/,a)', ' <<<+- DAMASK_spectral init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- DAMASK_grid init -+>>>'; flush(IO_STDOUT)
print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
!--------------------------------------------------------------------------------------------------
! initialize field solver information
nActiveFields = 1
if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1
if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1
allocate(solres(nActiveFields))
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! reading field paramters from numerics file and do sanity checks ! reading field paramters from numerics file and do sanity checks
@ -133,19 +131,22 @@ program DAMASK_grid
if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter')
if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack')
config_load => YAML_parse_file(trim(interface_loadFile))
solver => config_load%get('solver')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assign mechanics solver depending on selected type ! assign mechanics solver depending on selected type
debug_grid => config_debug%get('grid',defaultVal=emptyList) nActiveFields = 1
select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) select case (solver%get_asString('mechanical'))
case ('Basic') case ('spectral_basic')
mechanical_init => grid_mechanical_spectral_basic_init mechanical_init => grid_mechanical_spectral_basic_init
mechanical_forward => grid_mechanical_spectral_basic_forward mechanical_forward => grid_mechanical_spectral_basic_forward
mechanical_solution => grid_mechanical_spectral_basic_solution mechanical_solution => grid_mechanical_spectral_basic_solution
mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords
mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite
case ('Polarisation') case ('spectral_polarization')
mechanical_init => grid_mechanical_spectral_polarisation_init mechanical_init => grid_mechanical_spectral_polarisation_init
mechanical_forward => grid_mechanical_spectral_polarisation_forward mechanical_forward => grid_mechanical_spectral_polarisation_forward
mechanical_solution => grid_mechanical_spectral_polarisation_solution mechanical_solution => grid_mechanical_spectral_polarisation_solution
@ -160,32 +161,36 @@ program DAMASK_grid
mechanical_restartWrite => grid_mechanical_FEM_restartWrite mechanical_restartWrite => grid_mechanical_FEM_restartWrite
case default case default
call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver'))) call IO_error(error_ID = 891, ext_msg = trim(solver%get_asString('mechanical')))
end select end select
!--------------------------------------------------------------------------------------------------
! initialize field solver information
if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1
if (solver%get_asString('damage', defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1
allocate(solres(nActiveFields))
allocate( ID(nActiveFields))
field = 1
ID(field) = FIELD_MECH_ID ! mechanical active by default
thermalActive: if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') then
field = field + 1
ID(field) = FIELD_THERMAL_ID
endif thermalActive
damageActive: if (solver%get_asString('damage',defaultVal = 'n/a') == 'spectral') then
field = field + 1
ID(field) = FIELD_DAMAGE_ID
endif damageActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading information from load case file and to sanity checks
config_load => YAML_parse_file(trim(interface_loadFile))
load_steps => config_load%get('loadstep') load_steps => config_load%get('loadstep')
allocate(loadCases(load_steps%length)) ! array of load cases allocate(loadCases(load_steps%length)) ! array of load cases
do l = 1, load_steps%length do l = 1, load_steps%length
allocate(loadCases(l)%ID(nActiveFields))
field = 1
loadCases(l)%ID(field) = FIELD_MECH_ID ! mechanical active by default
thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
field = field + 1
loadCases(l)%ID(field) = FIELD_THERMAL_ID
endif thermalActive
damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then
field = field + 1
loadCases(l)%ID(field) = FIELD_DAMAGE_ID
endif damageActive
load_step => load_steps%get(l) load_step => load_steps%get(l)
step_bc => load_step%get('boundary_conditions') step_bc => load_step%get('boundary_conditions')
step_mech => step_bc%get('mechanical') step_mech => step_bc%get('mechanical')
@ -220,19 +225,17 @@ program DAMASK_grid
if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/dot_F/F missing') if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/dot_F/F missing')
step_discretization => load_step%get('discretization') step_discretization => load_step%get('discretization')
if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing') if (.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing')
if(.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing') if (.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing')
loadCases(l)%t = step_discretization%get_asFloat('t') loadCases(l)%t = step_discretization%get_asFloat('t')
loadCases(l)%N = step_discretization%get_asInt ('N') loadCases(l)%N = step_discretization%get_asInt ('N')
loadCases(l)%r = step_discretization%get_asFloat('r', defaultVal= 1.0_pReal) 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)) loadCases(l)%f_restart = load_step%get_asInt('f_restart', defaultVal=huge(0))
loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1) loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1)
loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. & loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1)
merge(.true.,.false.,l > 1))
reportAndCheck: if (worldrank == 0) then reportAndCheck: if (worldrank == 0) then
write (loadcase_string, '(i0)' ) l
print'(/,a,i0)', ' load case: ', l print'(/,a,i0)', ' load case: ', l
print*, ' estimate_rate:', loadCases(l)%estimate_rate print*, ' estimate_rate:', loadCases(l)%estimate_rate
if (loadCases(l)%deformation%myType == 'L') then if (loadCases(l)%deformation%myType == 'L') then
@ -292,7 +295,7 @@ program DAMASK_grid
if (loadCases(l)%f_restart < huge(0)) & if (loadCases(l)%f_restart < huge(0)) &
print'(a,i0)', ' f_restart: ', loadCases(l)%f_restart print'(a,i0)', ' f_restart: ', loadCases(l)%f_restart
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message if (errorID > 0) call IO_error(error_ID = errorID, el = l)
endif reportAndCheck endif reportAndCheck
enddo enddo
@ -301,12 +304,14 @@ program DAMASK_grid
! doing initialization depending on active solvers ! doing initialization depending on active solvers
call spectral_Utilities_init call spectral_Utilities_init
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(1)%ID(field)) select case (ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call mechanical_init call mechanical_init
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
call grid_thermal_spectral_init initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict)
thermal => initial_conditions%get('thermal',defaultVal=emptyDict)
call grid_thermal_spectral_init(thermal%get_asFloat('T',defaultVal = T_0))
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
call grid_damage_spectral_init call grid_damage_spectral_init
@ -377,7 +382,7 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forward fields ! forward fields
do field = 1, nActiveFields do field = 1, nActiveFields
select case(loadCases(l)%ID(field)) select case(ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call mechanical_forward (& call mechanical_forward (&
cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
@ -397,7 +402,7 @@ program DAMASK_grid
stagIterate = .true. stagIterate = .true.
do while (stagIterate) do while (stagIterate)
do field = 1, nActiveFields do field = 1, nActiveFields
select case(loadCases(l)%ID(field)) select case(ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
solres(field) = mechanical_solution(incInfo) solres(field) = mechanical_solution(incInfo)
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)

2
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -116,7 +116,7 @@ subroutine grid_mechanical_spectral_polarisation_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mechanical_spectral_polarisation init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mechanical_spectral_polarization init -+>>>'; flush(IO_STDOUT)
print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'

7
src/grid/grid_thermal_spectral.f90
View File

@ -61,7 +61,9 @@ contains
!> @brief allocates all neccessary fields and fills them with data !> @brief allocates all neccessary fields and fills them with data
! ToDo: Restart not implemented ! ToDo: Restart not implemented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_thermal_spectral_init subroutine grid_thermal_spectral_init(T_0)
real(pReal), intent(in) :: T_0
PetscInt, dimension(0:worldsize-1) :: localK PetscInt, dimension(0:worldsize-1) :: localK
integer :: i, j, k, ce integer :: i, j, k, ce
@ -131,9 +133,10 @@ subroutine grid_thermal_spectral_init
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
T_current(i,j,k) = homogenization_thermal_T(ce) T_current(i,j,k) = T_0
T_lastInc(i,j,k) = T_current(i,j,k) T_lastInc(i,j,k) = T_current(i,j,k)
T_stagInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k)
call homogenization_thermal_setField(T_0,0.0_pReal,ce)
enddo; enddo; enddo enddo; enddo; enddo
call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current

49
src/homogenization.f90
View File

@ -35,15 +35,11 @@ module homogenization
homogState, & homogState, &
damageState_h damageState_h
integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable :: &
real(pReal), dimension(:), allocatable, public, protected :: &
thermal_initialT
integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable, public, protected :: &
thermal_type !< thermal transport model thermal_type !< thermal transport model
integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: &
damage_type !< nonlocal damage model damage_type !< nonlocal damage model
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: &
homogenization_type !< type of each homogenization homogenization_type !< type of each homogenization
type, private :: tNumerics_damage type, private :: tNumerics_damage
@ -133,10 +129,8 @@ module homogenization
module function thermal_conduction_getConductivity(ce) result(K) module function thermal_conduction_getConductivity(ce) result(K)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), dimension(3,3) :: K real(pReal), dimension(3,3) :: K
end function thermal_conduction_getConductivity end function thermal_conduction_getConductivity
module function thermal_conduction_getSpecificHeat(ce) result(c_P) module function thermal_conduction_getSpecificHeat(ce) result(c_P)
@ -177,7 +171,6 @@ module homogenization
end function damage_nonlocal_getMobility end function damage_nonlocal_getMobility
module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
@ -185,21 +178,17 @@ module homogenization
phiDot, dPhiDot_dPhi phiDot, dPhiDot_dPhi
end subroutine damage_nonlocal_getSourceAndItsTangent end subroutine damage_nonlocal_getSourceAndItsTangent
module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
end subroutine damage_nonlocal_putNonLocalDamage end subroutine damage_nonlocal_putNonLocalDamage
module subroutine damage_nonlocal_results(ho,group) module subroutine damage_nonlocal_results(ho,group)
integer, intent(in) :: ho integer, intent(in) :: ho
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
end subroutine damage_nonlocal_results end subroutine damage_nonlocal_results
end interface end interface
public :: & public :: &
@ -242,21 +231,18 @@ subroutine homogenization_init()
allocate(homogState (size(material_name_homogenization))) allocate(homogState (size(material_name_homogenization)))
allocate(damageState_h (size(material_name_homogenization))) allocate(damageState_h (size(material_name_homogenization)))
call material_parseHomogenization call material_parseHomogenization()
num_homog => config_numerics%get('homogenization',defaultVal=emptyDict) num_homog => config_numerics%get('homogenization',defaultVal=emptyDict)
num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10) num%nMPstate = num_homogGeneric%get_asInt('nMPstate',defaultVal=10)
if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate') if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate')
call mechanical_init(num_homog) call mechanical_init(num_homog)
call thermal_init() call thermal_init()
call damage_init() call damage_init()
call damage_nonlocal_init()
call damage_nonlocal_init
end subroutine homogenization_init end subroutine homogenization_init
@ -326,7 +312,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
enddo enddo
!$OMP END DO !$OMP END DO
if (.not. terminallyIll ) then if (.not. terminallyIll) then
!$OMP DO PRIVATE(ho,ph,ce) !$OMP DO PRIVATE(ho,ph,ce)
do el = FEsolving_execElem(1),FEsolving_execElem(2) do el = FEsolving_execElem(1),FEsolving_execElem(2)
if (terminallyIll) continue if (terminallyIll) continue
@ -341,8 +327,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
print*, ' Integration point ', ip,' at element ', el, ' terminally ill' print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
terminallyIll = .true. ! ...and kills all others terminallyIll = .true. ! ...and kills all others
endif endif
call thermal_homogenize(ip,el)
enddo enddo
call thermal_homogenize(ip,el)
enddo enddo
enddo enddo
!$OMP END DO !$OMP END DO
@ -556,7 +542,6 @@ subroutine material_parseHomogenization
allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID)
allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID)
allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal)
do h=1, size(material_name_homogenization) do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h) homog => material_homogenization%get(h)
@ -572,27 +557,20 @@ subroutine material_parseHomogenization
call IO_error(500,ext_msg=homogMech%get_asString('type')) call IO_error(500,ext_msg=homogMech%get_asString('type'))
end select end select
if (homog%contains('thermal')) then
if(homog%contains('thermal')) then
homogThermal => homog%get('thermal') homogThermal => homog%get('thermal')
thermal_initialT(h) = homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal)
select case (homogThermal%get_asString('type')) select case (homogThermal%get_asString('type'))
case('isothermal') case('pass')
thermal_type(h) = THERMAL_isothermal_ID
case('conduction')
thermal_type(h) = THERMAL_conduction_ID thermal_type(h) = THERMAL_conduction_ID
case default case default
call IO_error(500,ext_msg=homogThermal%get_asString('type')) call IO_error(500,ext_msg=homogThermal%get_asString('type'))
end select end select
endif endif
if(homog%contains('damage')) then if (homog%contains('damage')) then
homogDamage => homog%get('damage') homogDamage => homog%get('damage')
select case (homogDamage%get_asString('type')) select case (homogDamage%get_asString('type'))
case('none') case('pass')
damage_type(h) = DAMAGE_none_ID
case('nonlocal')
damage_type(h) = DAMAGE_nonlocal_ID damage_type(h) = DAMAGE_nonlocal_ID
case default case default
call IO_error(500,ext_msg=homogDamage%get_asString('type')) call IO_error(500,ext_msg=homogDamage%get_asString('type'))
@ -600,7 +578,6 @@ subroutine material_parseHomogenization
endif endif
enddo enddo
end subroutine material_parseHomogenization end subroutine material_parseHomogenization

2
src/homogenization_thermal.f90
View File

@ -44,7 +44,7 @@ module subroutine thermal_init()
allocate(current(configHomogenizations%length)) allocate(current(configHomogenizations%length))
do ho = 1, configHomogenizations%length do ho = 1, configHomogenizations%length
allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=thermal_initialT(ho)) allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=300.0_pReal)
allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal) allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal)
configHomogenization => configHomogenizations%get(ho) configHomogenization => configHomogenizations%get(ho)
associate(prm => param(ho)) associate(prm => param(ho))

2
src/material.f90
View File

@ -30,7 +30,7 @@ module material
material_homogenizationAt, & !< homogenization ID of each element material_homogenizationAt, & !< homogenization ID of each element
material_homogenizationAt2, & !< per cell material_homogenizationAt2, & !< per cell
material_homogenizationMemberAt2 !< cell material_homogenizationMemberAt2 !< cell
integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem) integer, dimension(:,:), allocatable :: & ! (ip,elem)
material_homogenizationMemberAt !< position of the element within its homogenization instance material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt, & !< phase ID of each element material_phaseAt, & !< phase ID of each element

8
src/math.f90
View File

@ -1133,6 +1133,7 @@ real(pReal) pure function math_areaTriangle(v1,v2,v3)
real(pReal), dimension (3), intent(in) :: v1,v2,v3 real(pReal), dimension (3), intent(in) :: v1,v2,v3
math_areaTriangle = 0.5_pReal * norm2(math_cross(v1-v2,v1-v3)) math_areaTriangle = 0.5_pReal * norm2(math_cross(v1-v2,v1-v3))
end function math_areaTriangle end function math_areaTriangle
@ -1147,11 +1148,13 @@ real(pReal) pure elemental function math_clip(a, left, right)
real(pReal), intent(in) :: a real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: left, right real(pReal), intent(in), optional :: left, right
math_clip = a math_clip = a
if (present(left)) math_clip = max(left,math_clip) if (present(left)) math_clip = max(left,math_clip)
if (present(right)) math_clip = min(right,math_clip) if (present(right)) math_clip = min(right,math_clip)
if (present(left) .and. present(right)) & if (present(left) .and. present(right)) then
math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right) if(left>right) error stop 'left > right'
endif
end function math_clip end function math_clip
@ -1182,6 +1185,7 @@ subroutine selfTest
integer :: d integer :: d
logical :: e logical :: e
if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - &
math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) & math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) &
error stop 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]' error stop 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]'

8
src/phase_damage.f90
View File

@ -133,7 +133,7 @@ module subroutine damage_init
integer :: & integer :: &
ph, & !< counter in phase loop ph, & !< counter in phase loop
Nconstituents Nmembers
class(tNode), pointer :: & class(tNode), pointer :: &
phases, & phases, &
phase, & phase, &
@ -151,10 +151,10 @@ module subroutine damage_init
do ph = 1,phases%length do ph = 1,phases%length
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
allocate(current(ph)%phi(Nconstituents),source=1.0_pReal) allocate(current(ph)%phi(Nmembers),source=1.0_pReal)
allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal) allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal)
phase => phases%get(ph) phase => phases%get(ph)
sources => phase%get('damage',defaultVal=emptyList) sources => phase%get('damage',defaultVal=emptyList)

6
src/phase_damage_anisobrittle.f90
View File

@ -40,7 +40,7 @@ module function anisobrittle_init() result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Nconstituents,p integer :: Nmembers,p
integer, dimension(:), allocatable :: N_cl integer, dimension(:), allocatable :: N_cl
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -92,8 +92,8 @@ module function anisobrittle_init() result(mySources)
if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
Nconstituents = count(material_phaseAt==p) * discretization_nIPs Nmembers = count(material_phaseAt==p) * discretization_nIPs
call phase_allocateState(damageState(p),Nconstituents,1,1,0) call phase_allocateState(damageState(p),Nmembers,1,1,0)
damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'

6
src/phase_damage_anisoductile.f90
View File

@ -35,7 +35,7 @@ module function anisoductile_init() result(mySources)
pl, & pl, &
sources, & sources, &
src src
integer :: Ninstances,Nconstituents,p integer :: Ninstances,Nmembers,p
integer, dimension(:), allocatable :: N_sl integer, dimension(:), allocatable :: N_sl
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -78,8 +78,8 @@ module function anisoductile_init() result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
Nconstituents=count(material_phaseAt2==p) Nmembers=count(material_phaseAt2==p)
call phase_allocateState(damageState(p),Nconstituents,1,1,0) call phase_allocateState(damageState(p),Nmembers,1,1,0)
damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'

6
src/phase_damage_isobrittle.f90
View File

@ -31,7 +31,7 @@ module function isobrittle_init() result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Nconstituents,p integer :: Nmembers,p
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -64,8 +64,8 @@ module function isobrittle_init() result(mySources)
! sanity checks ! sanity checks
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
Nconstituents = count(material_phaseAt2==p) Nmembers = count(material_phaseAt2==p)
call phase_allocateState(damageState(p),Nconstituents,1,1,1) call phase_allocateState(damageState(p),Nmembers,1,1,1)
damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'

6
src/phase_damage_isoductile.f90
View File

@ -33,7 +33,7 @@ module function isoductile_init() result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstances,Nconstituents,p integer :: Ninstances,Nmembers,p
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -68,8 +68,8 @@ module function isoductile_init() result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
Nconstituents=count(material_phaseAt2==p) Nmembers=count(material_phaseAt2==p)
call phase_allocateState(damageState(p),Nconstituents,1,1,0) call phase_allocateState(damageState(p),Nmembers,1,1,0)
damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'

34
src/phase_mechanical.f90
View File

@ -191,7 +191,7 @@ module subroutine mechanical_init(materials,phases)
ph, & ph, &
me, & me, &
stiffDegradationCtr, & stiffDegradationCtr, &
Nconstituents Nmembers
class(tNode), pointer :: & class(tNode), pointer :: &
num_crystallite, & num_crystallite, &
material, & material, &
@ -229,22 +229,22 @@ module subroutine mechanical_init(materials,phases)
allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseMemberAt))) allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseMemberAt)))
do ph = 1, phases%length do ph = 1, phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs Nmembers = count(material_phaseAt == ph) * discretization_nIPs
allocate(phase_mechanical_Fi(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fe(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fi0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Lp0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Lp(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_F(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_F0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
phase => phases%get(ph) phase => phases%get(ph)
mech => phase%get('mechanics') mech => phase%get('mechanics')
@ -279,7 +279,6 @@ module subroutine mechanical_init(materials,phases)
endif endif
!$OMP PARALLEL DO PRIVATE(ph,me,material,constituents,constituent)
do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2)
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
material => materials%get(discretization_materialAt(el)) material => materials%get(discretization_materialAt(el))
@ -305,7 +304,6 @@ module subroutine mechanical_init(materials,phases)
enddo enddo
enddo; enddo enddo; enddo
!$OMP END PARALLEL DO
! initialize plasticity ! initialize plasticity

14
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -79,7 +79,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
ph, i, & ph, i, &
Nconstituents, & Nmembers, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -220,18 +220,18 @@ module function plastic_dislotungsten_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:) stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob = spread(rho_mob_0,2,Nconstituents) stt%rho_mob = spread(rho_mob_0,2,Nmembers)
dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:) dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) 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' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@ -239,7 +239,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:) stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip = spread(rho_dip_0,2,Nconstituents) stt%rho_dip = spread(rho_dip_0,2,Nmembers)
dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:) dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
@ -251,8 +251,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
! global alias ! global alias
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
allocate(dst%Lambda_sl(prm%sum_N_sl,Nconstituents), source=0.0_pReal) allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal)
allocate(dst%threshold_stress(prm%sum_N_sl,Nconstituents), source=0.0_pReal) allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal)
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally

143
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -24,7 +24,6 @@ submodule(phase:plastic) dislotwin
q_sb = 1.0_pReal, & !< q-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity
D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance
i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning
tau_0 = 1.0_pReal, & !< strength due to elements in solid solution
L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors
L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors
x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus
@ -53,6 +52,7 @@ submodule(phase:plastic) dislotwin
q, & !< q-exponent in glide velocity q, & !< q-exponent in glide velocity
r, & !< r-exponent in twin nucleation rate r, & !< r-exponent in twin nucleation rate
s, & !< s-exponent in trans nucleation rate s, & !< s-exponent in trans nucleation rate
tau_0, & !< strength due to elements in solid solution
gamma_char, & !< characteristic shear for twins gamma_char, & !< characteristic shear for twins
B !< drag coefficient B !< drag coefficient
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
@ -81,7 +81,7 @@ submodule(phase:plastic) dislotwin
logical :: & logical :: &
ExtendedDislocations, & !< consider split into partials for climb calculation ExtendedDislocations, & !< consider split into partials for climb calculation
fccTwinTransNucleation, & !< twinning and transformation models are for fcc fccTwinTransNucleation, & !< twinning and transformation models are for fcc
dipoleFormation !< flag indicating consideration of dipole formation omitDipoles !< flag controlling consideration of dipole formation
end type !< container type for internal constitutive parameters end type !< container type for internal constitutive parameters
type :: tDislotwinState type :: tDislotwinState
@ -127,7 +127,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
ph, i, & ph, i, &
Nconstituents, & Nmembers, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -213,10 +213,10 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl))
prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl)) prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl)) prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl))
prm%tau_0 = pl%get_asFloats('tau_0', requiredSize=size(N_sl))
prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), & prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal, i=1,size(N_sl))]) defaultVal=[(0.0_pReal, i=1,size(N_sl))])
prm%tau_0 = pl%get_asFloat('tau_0')
prm%D_a = pl%get_asFloat('D_a') prm%D_a = pl%get_asFloat('D_a')
prm%D_0 = pl%get_asFloat('D_0') prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl') prm%Q_cl = pl%get_asFloat('Q_cl')
@ -226,7 +226,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT') prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT')
endif endif
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation',defaultVal = .false.) prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.)
! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex)
! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 ! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
@ -242,6 +242,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%i_sl = math_expand(prm%i_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%p = math_expand(prm%p, N_sl) prm%p = math_expand(prm%p, N_sl)
prm%q = math_expand(prm%q, N_sl) prm%q = math_expand(prm%q, N_sl)
prm%tau_0 = math_expand(prm%tau_0, N_sl)
prm%B = math_expand(prm%B, N_sl) prm%B = math_expand(prm%B, N_sl)
! sanity checks ! sanity checks
@ -405,21 +406,21 @@ module function plastic_dislotwin_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
+ size(['f_tw']) * prm%sum_N_tw & + size(['f_tw']) * prm%sum_N_tw &
+ size(['f_tr']) * prm%sum_N_tr + size(['f_tr']) * prm%sum_N_tr
sizeState = sizeDotState sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! locally defined state aliases and initialization of state0 and atol
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:) stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob= spread(rho_mob_0,2,Nconstituents) stt%rho_mob= spread(rho_mob_0,2,Nmembers)
dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:) dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) 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' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@ -427,7 +428,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:) stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip= spread(rho_dip_0,2,Nconstituents) stt%rho_dip= spread(rho_dip_0,2,Nmembers)
dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:) dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
@ -443,28 +444,28 @@ module function plastic_dislotwin_init() result(myPlasticity)
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:) stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:) dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tw',defaultVal=1.0e-7_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tw'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tr endIndex = endIndex + prm%sum_N_tr
stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:) stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:) dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tr',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tr'
allocate(dst%Lambda_sl (prm%sum_N_sl,Nconstituents),source=0.0_pReal) allocate(dst%Lambda_sl (prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%tau_pass (prm%sum_N_sl,Nconstituents),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%Lambda_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal) allocate(dst%Lambda_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal)
allocate(dst%tau_hat_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal) allocate(dst%tau_hat_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal)
allocate(dst%tau_r_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal) allocate(dst%tau_r_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal)
allocate(dst%V_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal) allocate(dst%V_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal)
allocate(dst%Lambda_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) allocate(dst%Lambda_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal)
allocate(dst%tau_hat_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) allocate(dst%tau_hat_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal)
allocate(dst%tau_r_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) allocate(dst%tau_r_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal)
allocate(dst%V_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) allocate(dst%V_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal)
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
@ -535,9 +536,9 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_sl,ddot_gamma_dtau_slip dot_gamma_sl,ddot_gamma_dtau_slip
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_twin,ddot_gamma_dtau_twin dot_gamma_tw,ddot_gamma_dtau_tw
real(pReal), dimension(param(ph)%sum_N_tr) :: & real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr,ddot_gamma_dtau_trans dot_gamma_tr,ddot_gamma_dtau_tr
real(pReal):: dot_gamma_sb real(pReal):: dot_gamma_sb
real(pReal), dimension(3,3) :: eigVectors, P_sb real(pReal), dimension(3,3) :: eigVectors, P_sb
real(pReal), dimension(3) :: eigValues real(pReal), dimension(3) :: eigValues
@ -578,20 +579,20 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+ ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
enddo slipContribution enddo slipContribution
call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin,ddot_gamma_dtau_twin) call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw,ddot_gamma_dtau_tw)
twinContibution: do i = 1, prm%sum_N_tw twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinContibution enddo twinContibution
call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_trans) call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_tr)
transContibution: do i = 1, prm%sum_N_tr transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
enddo transContibution enddo transContibution
Lp = Lp * f_unrotated Lp = Lp * f_unrotated
@ -646,7 +647,6 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
f_unrotated, & f_unrotated, &
rho_dip_distance, & rho_dip_distance, &
v_cl, & !< climb velocity v_cl, & !< climb velocity
Gamma, & !< stacking fault energy
tau, & tau, &
sigma_cl, & !< climb stress sigma_cl, & !< climb stress
b_d !< ratio of Burgers vector to stacking fault width b_d !< ratio of Burgers vector to stacking fault width
@ -656,7 +656,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
rho_dip_distance_min, & rho_dip_distance_min, &
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_twin dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tr) :: & real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr dot_gamma_tr
@ -675,7 +675,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
slipState: do i = 1, prm%sum_N_sl slipState: do i = 1, prm%sum_N_sl
tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
significantSlipStress: if (dEq0(tau)) then significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then
dot_rho_dip_formation(i) = 0.0_pReal dot_rho_dip_formation(i) = 0.0_pReal
dot_rho_dip_climb(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal
else significantSlipStress else significantSlipStress
@ -683,24 +683,18 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me)) rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me))
rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i)) rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i))
if (prm%dipoleFormation) then
dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) & dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) &
* stt%rho_mob(i,me)*abs(dot_gamma_sl(i)) * stt%rho_mob(i,me)*abs(dot_gamma_sl(i))
else
dot_rho_dip_formation(i) = 0.0_pReal
endif
if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then
dot_rho_dip_climb(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal
else else
!@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
if (prm%ExtendedDislocations) then b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) &
Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T * (prm%Gamma_sf_0K + prm%dGamma_sf_dT * T) / (prm%mu*prm%b_sl(i)), &
b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) 1.0_pReal, &
else prm%ExtendedDislocations)
b_d = 1.0_pReal
endif
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & 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) * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
@ -718,8 +712,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) & - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) &
- dot_rho_dip_climb - dot_rho_dip_climb
call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin) call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw)
dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char dot%f_tw(:,me) = f_unrotated*dot_gamma_tw/prm%gamma_char
call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr) call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr)
dot%f_tr(:,me) = f_unrotated*dot_gamma_tr dot%f_tr(:,me) = f_unrotated*dot_gamma_tr
@ -741,11 +735,9 @@ module subroutine dislotwin_dependentState(T,ph,me)
T T
real(pReal) :: & real(pReal) :: &
sumf_twin,Gamma,sumf_trans sumf_tw,Gamma,sumf_tr
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation inv_lambda_sl
inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
f_over_t_tw f_over_t_tw
@ -760,38 +752,27 @@ module subroutine dislotwin_dependentState(T,ph,me)
stt => state(ph),& stt => state(ph),&
dst => dependentState(ph)) dst => dependentState(ph))
sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,me)) sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,me))
sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,me)) sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,me))
Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
!* rescaled volume fraction for topology !* rescaled volume fraction for topology
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw ! this is per system ... f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw ! this is per system ...
f_over_t_tr = sumf_trans/prm%t_tr ! but this not f_over_t_tr = sumf_tr/prm%t_tr ! but this not
! ToDo ...Physically correct, but naming could be adjusted ! ToDo ...Physically correct, but naming could be adjusted
inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, & inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl
stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl
if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
dst%Lambda_sl(:,me) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl)
inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans)
if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: better logic needed here
dst%Lambda_sl(:,me) = prm%D &
/ (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr))
else
dst%Lambda_sl(:,me) = prm%D &
/ (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct?
endif
inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
!* threshold stress for dislocation motion !* threshold stress for dislocation motion
@ -957,7 +938,7 @@ end subroutine kinetics_slip
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
dot_gamma_twin,ddot_gamma_dtau_twin) dot_gamma_tw,ddot_gamma_dtau_tw)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -970,9 +951,9 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
dot_gamma_twin dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: &
ddot_gamma_dtau_twin ddot_gamma_dtau_tw
real, dimension(param(ph)%sum_N_tw) :: & real, dimension(param(ph)%sum_N_tw) :: &
tau, & tau, &
@ -1004,16 +985,16 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
significantStress: where(tau > tol_math_check) significantStress: where(tau > tol_math_check)
StressRatio_r = (dst%tau_hat_tw(:,me)/tau)**prm%r StressRatio_r = (dst%tau_hat_tw(:,me)/tau)**prm%r
dot_gamma_twin = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r) dot_gamma_tw = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r)
ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
else where significantStress else where significantStress
dot_gamma_twin = 0.0_pReal dot_gamma_tw = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal ddot_gamma_dtau = 0.0_pReal
end where significantStress end where significantStress
end associate end associate
if(present(ddot_gamma_dtau_twin)) ddot_gamma_dtau_twin = ddot_gamma_dtau if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau
end subroutine kinetics_twin end subroutine kinetics_twin
@ -1026,7 +1007,7 @@ end subroutine kinetics_twin
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
dot_gamma_tr,ddot_gamma_dtau_trans) dot_gamma_tr,ddot_gamma_dtau_tr)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
@ -1041,7 +1022,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: &
dot_gamma_tr dot_gamma_tr
real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: &
ddot_gamma_dtau_trans ddot_gamma_dtau_tr
real, dimension(param(ph)%sum_N_tr) :: & real, dimension(param(ph)%sum_N_tr) :: &
tau, & tau, &
@ -1081,7 +1062,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
end associate end associate
if(present(ddot_gamma_dtau_trans)) ddot_gamma_dtau_trans = ddot_gamma_dtau if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau
end subroutine kinetics_trans end subroutine kinetics_trans

6
src/phase_mechanical_plastic_isotropic.f90
View File

@ -52,7 +52,7 @@ module function plastic_isotropic_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
ph, & ph, &
Nconstituents, & Nmembers, &
sizeState, sizeDotState sizeState, sizeDotState
real(pReal) :: & real(pReal) :: &
xi_0 !< initial critical stress xi_0 !< initial critical stress
@ -119,11 +119,11 @@ module function plastic_isotropic_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
sizeDotState = size(['xi ','gamma']) sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

8
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -62,7 +62,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
ph, o, & ph, o, &
Nconstituents, & Nmembers, &
sizeState, sizeDeltaState, sizeDotState, & sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -165,19 +165,19 @@ module function plastic_kinehardening_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%crss => plasticState(ph)%state (startIndex:endIndex,:) stt%crss => plasticState(ph)%state (startIndex:endIndex,:)
stt%crss = spread(xi_0, 2, Nconstituents) stt%crss = spread(xi_0, 2, Nmembers)
dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:) dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) 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' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

61
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -177,7 +177,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
integer :: & integer :: &
Ninstances, & Ninstances, &
ph, & ph, &
Nconstituents, & Nmembers, &
sizeState, sizeDotState, sizeDependentState, sizeDeltaState, & sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
s1, s2, & s1, s2, &
s, t, l s, t, l
@ -398,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@ -412,9 +412,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState sizeDeltaState = sizeDotState
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
allocate(geom(ph)%V_0(Nconstituents)) allocate(geom(ph)%V_0(Nmembers))
call storeGeometry(ph) call storeGeometry(ph)
plasticState(ph)%nonlocal = pl%get_asBool('nonlocal') plasticState(ph)%nonlocal = pl%get_asBool('nonlocal')
@ -486,26 +486,26 @@ module function plastic_nonlocal_init() result(myPlasticity)
dot%rho_dip_scr => plasticState(ph)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) dot%rho_dip_scr => plasticState(ph)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho_dip_scr => plasticState(ph)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) del%rho_dip_scr => plasticState(ph)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%gamma => plasticState(ph)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) stt%gamma => plasticState(ph)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
dot%gamma => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) dot%gamma => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
del%gamma => plasticState(ph)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) del%gamma => plasticState(ph)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal) plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal)
if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) & if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma' extmsg = trim(extmsg)//' atol_gamma'
plasticState(ph)%slipRate => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) plasticState(ph)%slipRate => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents) stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers)
stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents) stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents) stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers)
stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents) stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers)
stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents) stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers)
stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents) stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
allocate(dst%tau_pass(prm%sum_N_sl,Nconstituents),source=0.0_pReal) allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal) allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal)
end associate end associate
if (Nconstituents > 0) call stateInit(ini,ph,Nconstituents) if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
plasticState(ph)%state0 = plasticState(ph)%state plasticState(ph)%state0 = plasticState(ph)%state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -527,7 +527,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
if(.not. myPlasticity(ph)) cycle if(.not. myPlasticity(ph)) cycle
phase => phases%get(ph) phase => phases%get(ph)
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
l = 0 l = 0
do t = 1,4 do t = 1,4
do s = 1,param(ph)%sum_N_sl do s = 1,param(ph)%sum_N_sl
@ -1579,13 +1579,13 @@ end subroutine plastic_nonlocal_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief populates the initial dislocation density !> @brief populates the initial dislocation density
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine stateInit(ini,phase,Nconstituents) subroutine stateInit(ini,phase,Nmembers)
type(tInitialParameters) :: & type(tInitialParameters) :: &
ini ini
integer,intent(in) :: & integer,intent(in) :: &
phase, & phase, &
Nconstituents Nmembers
integer :: & integer :: &
i, & i, &
e, & e, &
@ -1602,7 +1602,7 @@ subroutine stateInit(ini,phase,Nconstituents)
totalVolume, & totalVolume, &
densityBinning, & densityBinning, &
minimumIpVolume minimumIpVolume
real(pReal), dimension(Nconstituents) :: & real(pReal), dimension(Nmembers) :: &
volume volume
@ -1622,13 +1622,13 @@ subroutine stateInit(ini,phase,Nconstituents)
meanDensity = 0.0_pReal meanDensity = 0.0_pReal
do while(meanDensity < ini%random_rho_u) do while(meanDensity < ini%random_rho_u)
call random_number(rnd) call random_number(rnd)
phasemember = nint(rnd(1)*real(Nconstituents,pReal) + 0.5_pReal) phasemember = nint(rnd(1)*real(Nmembers,pReal) + 0.5_pReal)
s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume
stt%rhoSglMobile(s,phasemember) = densityBinning stt%rhoSglMobile(s,phasemember) = densityBinning
enddo enddo
else ! homogeneous distribution with noise else ! homogeneous distribution with noise
do e = 1, Nconstituents do e = 1, Nmembers
do f = 1,size(ini%N_sl,1) do f = 1,size(ini%N_sl,1)
from = 1 + sum(ini%N_sl(1:f-1)) from = 1 + sum(ini%N_sl(1:f-1))
upto = sum(ini%N_sl(1:f)) upto = sum(ini%N_sl(1:f))
@ -1809,7 +1809,7 @@ pure function getRho0(ph,me)
getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)
getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal)
where(abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & where (abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho0 = 0.0_pReal getRho0 = 0.0_pReal
end associate end associate
@ -1822,16 +1822,13 @@ subroutine storeGeometry(ph)
integer, intent(in) :: ph integer, intent(in) :: ph
integer :: ip, el, ce, co integer :: ip, el, ce, co
real(pReal), dimension(:), allocatable :: V
ce = 0
do el = 1, size(material_homogenizationMemberAt,2) V = reshape(IPvolume,[product(shape(IPvolume))])
do ip = 1, size(material_homogenizationMemberAt,1) do ce = 1, size(material_homogenizationMemberAt2,1)
ce = ce + 1
do co = 1, homogenization_maxNconstituents do co = 1, homogenization_maxNconstituents
if(material_phaseAt2(co,ce) == ph) then if (material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce)
geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = IPvolume(ip,el)
endif
enddo
enddo enddo
enddo enddo

10
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -71,7 +71,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
ph, i, & ph, i, &
Nconstituents, & Nmembers, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -223,20 +223,20 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(ph)%state (startIndex:endIndex,:) stt%xi_slip => plasticState(ph)%state (startIndex:endIndex,:)
stt%xi_slip = spread(xi_0_sl, 2, Nconstituents) stt%xi_slip = spread(xi_0_sl, 2, Nmembers)
dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:) dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) 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' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@ -244,7 +244,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(ph)%state (startIndex:endIndex,:) stt%xi_twin => plasticState(ph)%state (startIndex:endIndex,:)
stt%xi_twin = spread(xi_0_tw, 2, Nconstituents) stt%xi_twin = spread(xi_0_tw, 2, Nmembers)
dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:) dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) 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' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

38
src/phase_thermal.f90
View File

@ -15,13 +15,13 @@ submodule(phase) thermal
THERMAL_EXTERNALHEAT_ID THERMAL_EXTERNALHEAT_ID
end enum end enum
type :: tDataContainer type :: tDataContainer ! ?? not very telling name. Better: "fieldQuantities" ??
real(pReal), dimension(:), allocatable :: T, dot_T real(pReal), dimension(:), allocatable :: T, dot_T
end type tDataContainer end type tDataContainer
integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: &
thermal_source thermal_source
type(tDataContainer), dimension(:), allocatable :: current type(tDataContainer), dimension(:), allocatable :: current ! ?? not very telling name. Better: "field" ??
integer :: thermal_source_maxSizeDotState integer :: thermal_source_maxSizeDotState
@ -78,41 +78,36 @@ module subroutine thermal_init(phases)
integer :: & integer :: &
ph, so, & ph, so, &
Nconstituents Nmembers
print'(/,a)', ' <<<+- phase:thermal init -+>>>' print'(/,a)', ' <<<+- phase:thermal init -+>>>'
allocate(current(phases%length)) allocate(current(phases%length))
allocate(thermalState (phases%length)) allocate(thermalState(phases%length))
allocate(thermal_Nsources(phases%length),source = 0) allocate(thermal_Nsources(phases%length),source = 0)
do ph = 1, phases%length do ph = 1, phases%length
Nmembers = count(material_phaseAt2 == ph)
Nconstituents = count(material_phaseAt2 == ph) allocate(current(ph)%T(Nmembers),source=300.0_pReal)
allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
allocate(current(ph)%T(Nconstituents),source=300.0_pReal)
allocate(current(ph)%dot_T(Nconstituents),source=0.0_pReal)
phase => phases%get(ph) phase => phases%get(ph)
if(phase%contains('thermal')) then thermal => phase%get('thermal',defaultVal=emptyDict)
thermal => phase%get('thermal')
sources => thermal%get('source',defaultVal=emptyList) sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length thermal_Nsources(ph) = sources%length
endif
allocate(thermalstate(ph)%p(thermal_Nsources(ph))) allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
enddo enddo
allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
if(maxval(thermal_Nsources) /= 0) then if (maxval(thermal_Nsources) /= 0) then
where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
endif endif
thermal_source_maxSizeDotState = 0 thermal_source_maxSizeDotState = 0
PhaseLoop2:do ph = 1,phases%length do ph = 1,phases%length
do so = 1,thermal_Nsources(ph) do so = 1,thermal_Nsources(ph)
thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0 thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0
@ -120,7 +115,7 @@ module subroutine thermal_init(phases)
thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, &
maxval(thermalState(ph)%p%sizeDotState)) maxval(thermalState(ph)%p%sizeDotState))
enddo PhaseLoop2 enddo
end subroutine thermal_init end subroutine thermal_init
@ -156,7 +151,6 @@ module subroutine phase_thermal_getRate(TDot, ph,me)
Tdot = Tdot + my_Tdot Tdot = Tdot + my_Tdot
enddo enddo
end subroutine phase_thermal_getRate end subroutine phase_thermal_getRate
@ -185,7 +179,7 @@ function phase_thermal_collectDotState(ph,me) result(broken)
end function phase_thermal_collectDotState end function phase_thermal_collectDotState
module function thermal_stress(Delta_t,ph,me) result(converged_) module function thermal_stress(Delta_t,ph,me) result(converged_) ! ?? why is this called "stress" when it seems closer to "updateState" ??
real(pReal), intent(in) :: Delta_t real(pReal), intent(in) :: Delta_t
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
@ -212,7 +206,7 @@ function integrateThermalState(Delta_t, ph,me) result(broken)
sizeDotState sizeDotState
broken = phase_thermal_collectDotState(ph,me) broken = phase_thermal_collectDotState(ph,me)
if(broken) return if (broken) return
do so = 1, thermal_Nsources(ph) do so = 1, thermal_Nsources(ph)
sizeDotState = thermalState(ph)%p(so)%sizeDotState sizeDotState = thermalState(ph)%p(so)%sizeDotState
@ -301,14 +295,12 @@ function thermal_active(source_label,src_length) result(active_source)
allocate(active_source(src_length,phases%length), source = .false. ) allocate(active_source(src_length,phases%length), source = .false. )
do p = 1, phases%length do p = 1, phases%length
phase => phases%get(p) phase => phases%get(p)
if (phase%contains('thermal')) then thermal => phase%get('thermal',defaultVal=emptyDict)
thermal => phase%get('thermal',defaultVal=emptyList)
sources => thermal%get('source',defaultVal=emptyList) sources => thermal%get('source',defaultVal=emptyList)
do s = 1, sources%length do s = 1, sources%length
src => sources%get(s) src => sources%get(s)
if(src%get_asString('type') == source_label) active_source(s,p) = .true. active_source(s,p) = src%get_asString('type') == source_label
enddo enddo
endif
enddo enddo

6
src/phase_thermal_dissipation.f90
View File

@ -31,7 +31,7 @@ module function dissipation_init(source_length) result(mySources)
phase, & phase, &
sources, thermal, & sources, thermal, &
src src
integer :: so,Nconstituents,ph integer :: so,Nmembers,ph
mySources = thermal_active('dissipation',source_length) mySources = thermal_active('dissipation',source_length)
@ -54,8 +54,8 @@ module function dissipation_init(source_length) result(mySources)
src => sources%get(so) src => sources%get(so)
prm%kappa = src%get_asFloat('kappa') prm%kappa = src%get_asFloat('kappa')
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
call phase_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0)
end associate end associate
endif endif

6
src/phase_thermal_externalheat.f90
View File

@ -38,7 +38,7 @@ module function externalheat_init(source_length) result(mySources)
phase, & phase, &
sources, thermal, & sources, thermal, &
src src
integer :: so,Nconstituents,ph integer :: so,Nmembers,ph
mySources = thermal_active('externalheat',source_length) mySources = thermal_active('externalheat',source_length)
@ -67,8 +67,8 @@ module function externalheat_init(source_length) result(mySources)
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
Nconstituents = count(material_phaseAt2 == ph) Nmembers = count(material_phaseAt2 == ph)
call phase_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0)
end associate end associate
endif endif
enddo enddo