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 python-polishing

This commit is contained in:
Martin Diehl 2022-02-16 23:13:46 +01:00
parent 3ee98164fa 138b0e03f2
commit 9fc6469b13
93 changed files with 1469 additions and 962 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 5774122bf48d637704bb4afb10b87c34a4dbcaba
Subproject commit 0d639a9ba41db279b0d2825c8e8eddf0ccd91326

5
examples/config/phase/AISI304.yaml
View File

@ -1,6 +1,7 @@
references:
- H.M. Ledbetter
physica status solidi (a) 85(1):89-96, 1984
- H.M. Ledbetter,
physica status solidi (a) 85(1):89-96, 1984,
https://doi.org/10.1002/pssa.2210850111
lattice: cF
rho: 7937.0

1
examples/config/phase/Ag.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Silver
lattice: cF
rho: 10490.0

1
examples/config/phase/Al.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Aluminium
lattice: cF
rho: 2700.0

1
examples/config/phase/Au.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Gold
lattice: cF
rho: 19300.0

1
examples/config/phase/Cu.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Copper
lattice: cF
rho: 8960.0

1
examples/config/phase/Fe.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Iron
lattice: cI
rho: 7874.0

1
examples/config/phase/Mg.yaml
View File

@ -2,6 +2,7 @@ references:
- D. Tromans,
International Journal of Recent Research and Applied Studies 6(4):462-483, 2011,
https://www.arpapress.com/Volumes/Vol6Issue4/IJRRAS_6_4_14.pdf
lattice: hP
c/a: 1.62350
rho: 1740.0

1
examples/config/phase/Ni.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Nickel
lattice: cF
rho: 8908.0

9
examples/config/phase/Sn-beta.yaml Normal file
View File

@ -0,0 +1,9 @@
references:
- J.A. Rayne and B.S. Chandrasekhar,
Physical Review 120(5):1658-1663, 1960,
https://doi.org/10.1103/PhysRev.120.1658
- https://en.wikipedia.org/wiki/Tin
lattice: tI
c/a: 0.5458 # T=300K (c=31.83nm, a=5.832nm)
rho: 7265.0

1
examples/config/phase/Ti.yaml
View File

@ -1,6 +1,7 @@
references:
- https://www.totalmateria.com/page.aspx?ID=CheckArticle&site=ktn&NM=221
- https://en.wikipedia.org/wiki/Titanium
lattice: hP
c/a: 1.587
rho: 4506.0

1
examples/config/phase/W.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Tungsten
lattice: cI
rho: 19300.0

2
examples/config/phase/bct.yaml Normal file
View File

@ -0,0 +1,2 @@
lattice: tI
c/a: 0.55

6
examples/config/phase/damage/anisobrittle_cubic.yaml
View File

@ -1,11 +1,13 @@
type: anisobrittle
output: [f_phi]
N_cl: [3]
g_crit: [0.5e+7]
s_crit: [0.006666]
dot_o: 1.e-3
q: 20
output: [f_phi]
D_11: 1.0
mu: 0.001

2
examples/config/phase/damage/isobrittle_generic.yaml
View File

@ -1,7 +1,7 @@
type: isobrittle
W_crit: 1400000.0
output: [f_phi]
W_crit: 1400000.0
D_11: 1.0
mu: 0.001

2
examples/config/phase/hcp.yaml Normal file
View File

@ -0,0 +1,2 @@
lattice: hP
c/a: 1.6333

16
examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml
View File

@ -1,7 +1,13 @@
type: thermalexpansion
references:
- R.H. Bogaard et al.
Thermochimica Acta 218:373-393, 1993
https://doi.org/10.1016/0040-6031(93)80437-F
A_11: 15.0e-6
T_ref: 300.0
- R.H. Bogaard et al.,
Thermochimica Acta 218:373-393, 1993,
https://doi.org/10.1016/0040-6031(93)80437-F,
fit to Fig. 6 (T_min=100K, T_max=1400K)
A_11: 2.068e-08
A_11,T: 1.579e-09
A_11,T^2: 3.449e-13
T_ref: 293.15

6
examples/config/phase/mechanical/eigen/thermalexpansion_Al.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 23.1e-6
T_ref: 293.15

6
examples/config/phase/mechanical/eigen/thermalexpansion_Au.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 14.e-6
T_ref: 293.15

5
examples/config/phase/mechanical/eigen/thermalexpansion_C35E.yaml
View File

@ -1,8 +1,11 @@
type: thermalexpansion
references:
- https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg,
fitted from image description (Scilab code)
fit to image description (Scilab code)
A_11: 12.70371e-6
A_11,T: 7.54e-9
A_11,T^2: -1.0e-11
T_ref: 273.0

6
examples/config/phase/mechanical/eigen/thermalexpansion_Cu.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 17.e-6
T_ref: 293.15

6
examples/config/phase/mechanical/eigen/thermalexpansion_Fe.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 11.8e-6
T_ref: 293.15

17
examples/config/phase/mechanical/eigen/thermalexpansion_Sn-beta.yaml Normal file
View File

@ -0,0 +1,17 @@
type: thermalexpansion
references:
- V.T. Deshpande and D.B. Sirdeshmukh,
Acta Crystallographica 15:294-295, 1962,
https://doi.org/10.1107/S0365110X62000742,
fit to Tab. 2 (T_min=30ºC, T_max=210ºC)
A_11: 1.639e-05
A_11,T: 1.799e-08
A_11,T^2: 1.734e-10
A_33: 3.263e-05
A_33,T: 1.387e-08
A_33,T^2: 5.794e-10
T_ref: 293.15

6
examples/config/phase/mechanical/eigen/thermalexpansion_W.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 4.5e-6
T_ref: 293.15

7
examples/config/phase/mechanical/eigen/thermalexpansion_X20Cr13.yaml
View File

@ -1,7 +1,10 @@
type: thermalexpansion
references:
- https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg
fitted from image description (Scilab code)
- https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg,
fit to image description (Scilab code)
A_11: 11.365e-6
A_11,T: 5.0e-9
T_ref: 273.0

6
examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml
View File

@ -1,8 +1,10 @@
type: Hooke
references:
- H.M. Ledbetter
physica status solidi (a) 85(1):89-96, 1984
- H.M. Ledbetter,
physica status solidi (a) 85(1):89-96, 1984,
https://doi.org/10.1002/pssa.2210850111
C_11: 204.6e+9
C_12: 137.7e+9
C_44: 126.2e+9

5
examples/config/phase/mechanical/elastic/Hooke_Ag.yaml
View File

@ -1,4 +1,5 @@
type: Hooke
references:
- J.R. Neighbours and G.A. Alers,
Physical Review 111:707-712, 1958,
@ -7,8 +8,6 @@ references:
Journal of Applied Physics 37:3567-3572, 1966,
https://doi.org/10.1063/1.1708903
T_ref: 300
C_11: 122.9e+9
C_11,T: -313.5e+5
C_11,T^2: -107.3e+2
@ -20,3 +19,5 @@ C_12,T^2: -681.6e+1
C_44: 42.63e+9
C_44,T: -180.5e+5
C_44,T^2: -353.8e+1
T_ref: 300

27
examples/config/phase/mechanical/elastic/Hooke_Al.yaml
View File

@ -1,22 +1,25 @@
type: Hooke
references:
- G.N. Kamm and G.A. Alers,
Journal of Applied Physics 35:327-330, 1964,
https://doi.org/10.1063/1.1713309
https://doi.org/10.1063/1.1713309,
fit to Tab. I (T_min=100K, T_max=300K)
- D. Gerlich and E.S. Fisher,
Journal of Physics and Chemistry of Solids 30:1197-1205, 1969
https://doi.org/10.1016/0022-3697(69)90377-1
https://doi.org/10.1016/0022-3697(69)90377-1,
fit to Tab. 2 (T_min=293K, T_max=900K)
T_ref: 300
C_11: 106.9e+9
C_11,T: -3.685e+7
C_11,T^2: -1.016e+4
C_11: 106.1e+9
C_11,T: -359.3e+5
C_11,T^2: -152.7e+2
C_12: 60.55e+9
C_12,T: -8.307e+6
C_12,T^2: -4.353e+3
C_12: 57.83e+9
C_12,T: -781.6e+4
C_12,T^2: -551.3e+1
C_44: 28.37e+9
C_44,T: -1.418e+7
C_44,T^2: -3.253e+3
C_44: 24.31e+9
C_44,T: -142.9e+5
C_44,T^2: -404.6e+1
T_ref: 293.15

2
examples/config/phase/mechanical/elastic/Hooke_Au.yaml
View File

@ -1,9 +1,11 @@
type: Hooke
references:
- J.P. Hirth and J. Lothe,
Theory of Dislocations, 1982,
John Wiley & Sons,
page 837
C_11: 186.e+9
C_12: 157.e+9
C_44: 42.e+9

24
examples/config/phase/mechanical/elastic/Hooke_Cu.yaml
View File

@ -1,7 +1,21 @@
type: Hooke
references:
- https://www.mit.edu/~6.777/matprops/copper.htm,
fixed typo
C_11: 168.3e+9
C_12: 122.1e+9
C_44: 75.7e+9
- W.C. Overton, Jr. and J. Gaffney,
Physical Review 98(4):969-977, 1955,
https://doi.org/10.1103/PhysRev.98.969,
fit to Tab. I (T_min=100K, T_max=300K)
C_11: 168.6e+9
C_11,T: -3.779e+7
C_11,T^2: -2.536e+4
C_12: 121.5e+9
C_12,T: -1.632e+7
C_12,T^2: -1.116e+4
C_44: 75.59e+9
C_44,T: -2.912e+7
C_44,T^2: -1.669e+4
T_ref: 293.15

22
examples/config/phase/mechanical/elastic/Hooke_Fe.yaml
View File

@ -1,19 +1,21 @@
type: Hooke
references:
- D.J. Dever,
Journal of Applied Physics 43(8):3293-3301, 1972,
https://doi.org/10.1063/1.1661710
fit to Tab. II (T_min=25ºC, T_max=880ºC)
T_ref: 300
C_11: 232.1e+9
C_11,T: -4.678e+7
C_11,T^2: -5.762e+4
C_11: 231.7e+9
C_11,T: -47.6e+6
C_11,T^2: -54.4e+3
C_12: 135.9e+9
C_12,T: -1.695e+7
C_12,T^2: 1.555e+3
C_12: 135.8e+9
C_12,T: -12.9e+6
C_12,T^2: -7.3e+3
C_44: 117.0e+9
C_44,T: -2.047e+7
C_44,T^2: -2.814e+2
C_44: 116.8e+9
C_44,T: -19.4e+6
C_44,T^2: -2.5e+3
T_ref: 293.15

35
examples/config/phase/mechanical/elastic/Hooke_Mg.yaml
View File

@ -1,10 +1,29 @@
type: Hooke
references:
- D. Tromans,
International Journal of Recent Research and Applied Studies 6(4):462-483, 2011,
https://www.arpapress.com/Volumes/Vol6Issue4/IJRRAS_6_4_14.pdf
C_11: 59.3e+9
C_33: 61.5e+9
C_44: 16.4e+9
C_12: 25.7e+9
C_13: 21.4e+9
- L.J. Slutsky and C.W. Garland,
Physical Review 107(4):972-976, 1957,
https://doi.org/10.1103/PhysRev.107.972,
fit to Tab. I (T_min=100K, T_max=300K)
C_11: 59.50e+9
C_11,T: -1.930e+7
C_11,T^2: -1.215e+4
C_33: 61.72e+9
C_33,T: -2.175e+7
C_33,T^2: -5.624e+3
C_44: 16.46e+9
C_44,T: -1.006e+7
C_44,T^2: -7.692e+3
C_12: 25.62e+9
C_12,T: -2.216e+6
C_12,T^2: -4.138e+3
C_13: 21.46e+9
C_13,T: -1.921e+6
C_13,T^2: -4.283e+3
T_ref: 293.15

26
examples/config/phase/mechanical/elastic/Hooke_Ni.yaml
View File

@ -1,9 +1,21 @@
type: Hooke
references:
- J.P. Hirth and J. Lothe,
Theory of Dislocations, 1982,
John Wiley & Sons,
page 837
C_11: 246.5e+9
C_12: 147.3e+9
C_44: 124.7e+9
- G.A. Alers,
Journal of Physics and Chemistry of Solids 13(1-2):40-55, 1960,
https://doi.org/10.1016/0022-3697(60)90125-6,
fit to Tab. 2 (T_min=100K, T_max=700K)
C_11: 251.0e+9
C_11,T: -4.998e+7
C_11,T^2: -2.952e+4
C_12: 150.0e+9
C_12,T: -4.269e+6
C_12,T^2: -6.429e+3
C_44: 123.7e+9
C_44,T: -3.618e+7
C_44,T^2: -7.024e+3
T_ref: 293.15

2
examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml
View File

@ -1,8 +1,10 @@
type: Hooke
references:
- S.A. Kim and W.L. Johnson,
Materials Science & Engineering A 452-453:633-639, 2007,
https://doi.org/10.1016/j.msea.2006.11.147
C_11: 268.1e+9
C_12: 111.2e+9
C_44: 79.06e+9

32
examples/config/phase/mechanical/elastic/Hooke_Sn-beta.yaml Normal file
View File

@ -0,0 +1,32 @@
type: Hooke
references:
- J.A. Rayne and B.S. Chandrasekhar,
Physical Review 120(5):1658-1663, 1960,
https://doi.org/10.1103/PhysRev.120.1658,
fit to Fig. 2 (T_min=100K, T_max=300K) and Tab. IV (C_13, T_min=77K, T_max=300K)
C_11: 72.90e+9
C_11,T: -4.399e+7
C_11,T^2: -2.645e+4
C_12: 59.27e+9
C_12,T: 1.058e+7
C_12,T^2: 1.002e+4
C_13: 35.80e+9
C_13,T: -2.870e+6
C_33: 88.78e+9
C_33,T: -5.250e+7
C_33,T^2: 3.546e+3
C_44: 22.26e+9
C_44,T: -1.982e+7
C_44,T^2: -8.711e+3
C_66: 24.18e+9
C_66,T: -1.806e+7
C_66,T^2: -4.112e+3
T_ref: 293.15

2
examples/config/phase/mechanical/elastic/Hooke_TWIP-steel.yaml
View File

@ -1,4 +1,5 @@
type: Hooke
references:
- D. Music et al.,
Applied Physics Letters 99(19):191904, 2007,
@ -6,6 +7,7 @@ references:
- S.L. Wong et al.,
Acta Materialia 118:140-151, 2016,
https://doi.org/10.1016/j.actamat.2016.07.032
C_11: 175.0e+9
C_12: 115.0e+9
C_44: 135.0e+9

39
examples/config/phase/mechanical/elastic/Hooke_Ti.yaml
View File

@ -1,10 +1,33 @@
type: Hooke
references:
- L. Wang et al.,
Acta Materialia 132:598-610, 2017,
https://doi.org/10.1016/j.actamat.2017.05.015
C_11: 162.4e+9
C_33: 181.6e+9
C_44: 47.2e+9
C_12: 92.e+9
C_13: 69.e+9
- E.S. Fisher and C.J. Renken,
Physical Review 135(2A):A482-A494, 1964,
https://doi.org/10.1103/PhysRev.135.A482,
fit to Tab. IV (T_min=150K, T_max=250K)
- H. Ogi et al.,
Acta Materialia 52(7):2075-2080, 2004,
https://doi.org/10.1016/j.actamat.2004.01.002,
fit to Fig. 3 (T_min=300K, T_max=900K)
C_11: 162.6e+9
C_11,T: -6.150e+7
C_11,T^2: -5.557e+2
C_33: 183.3e+9
C_33,T: -1.655e+07
C_33,T^2: -1.022e+04
C_44: 45.80e+9
C_44,T: -2.936e+07
C_44,T^2: 7.120e+02
C_12: 89.97e+9
C_12,T: 2.776e+6
C_12,T^2: -2.389e+4
C_13: 69.53e+9
C_13,T: 1.057e+7
C_13,T^2: -2.966e+3
T_ref: 293.15

24
examples/config/phase/mechanical/elastic/Hooke_W.yaml
View File

@ -1,8 +1,20 @@
type: Hooke
references:
- D. Cereceda et al.,
International Journal of Plasticity 78:242-265, 2016,
https://doi.org/10.1016/j.ijplas.2015.09.002
C_11: 523.e+9
C_12: 202.e+9
C_44: 161.e+9
- F.H. Featherston and J.R. Nieghbours,
Physical Review 130(4):1324-1333,
https://doi.org/10.1103/PhysRev.130.1324,
fit to Tab. III (T_min=100K, T_max=300K)
C_11: 523.6e+9
C_11,T: -7.607e+7
C_11,T^2: -1.551e+5
C_12: 205.1e+9
C_12,T: -2.843e+6
C_44: 160.8e+9
C_44,T: -1.057e+7
C_44,T^2: 9.933e+3
T_ref: 293.15

2
examples/config/phase/mechanical/elastic/Hooke_vanishing-Poisson-ratio.yaml
View File

@ -1,8 +1,10 @@
type: Hooke
references:
- T. Maiti and P. Eisenlohr,
Scripta Materialia 145:37-40, 2018,
https://doi.org/10.1016/j.scriptamat.2017.09.047
C_11: 1.e+8
C_12: 1.e+6
C_44: 4.95e+7

4
examples/config/phase/mechanical/plastic/dislotungsten_W.yaml
View File

@ -1,4 +1,5 @@
type: dislotungsten
references:
- D. Cereceda et al.,
International Journal of Plasticity 78:242-265, 2016,
@ -6,8 +7,11 @@ references:
- R. Gröger et al.,
Acta Materialia 56(19):5412-5425, 2008,
https://doi.org/10.1016/j.actamat.2008.07.037
output: [Lambda_sl]
N_sl: [12]
b_sl: [2.72e-10]
rho_mob_0: [1.0e+9] # estimated from section 3.2
rho_dip_0: [1.0] # not given

4
examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml
View File

@ -1,4 +1,5 @@
type: dislotwin
references:
- K. Sedighiani et al.,
International Journal of Plasticity 134:102779, 2020,
@ -6,8 +7,11 @@ references:
- K. Sedighiani et al.,
Mechanics of Materials, 164:104117, 2022,
https://doi.org/10.1016/j.mechmat.2021.104117
output: [rho_dip, rho_mob]
N_sl: [12, 12]
b_sl: [2.49e-10, 2.49e-10]
rho_mob_0: [2.81e12, 2.8e+12]
rho_dip_0: [1.0, 1.0] # not given

3
examples/config/phase/mechanical/plastic/isotropic_free-surface.yaml
View File

@ -1,9 +1,12 @@
type: isotropic
references:
- T. Maiti and P. Eisenlohr,
Scripta Materialia 145:37-40, 2018,
https://doi.org/10.1016/j.scriptamat.2017.09.047
output: [xi]
dot_gamma_0: 0.001
n: 20.
xi_0: 0.3e+6

11
examples/config/phase/mechanical/plastic/nonlocal_Al.yaml
View File

@ -1,10 +1,13 @@
type: nonlocal
references:
C. Kords,
On the role of dislocation transport in the constitutive description of crystal plasticity,
RWTH Aachen 2013,
http://publications.rwth-aachen.de/record/229993/files/4862.pdf
- C. Kords,
On the role of dislocation transport in the constitutive description of crystal plasticity,
RWTH Aachen 2013,
http://publications.rwth-aachen.de/record/229993/files/4862.pdf
output: [rho_u_ed_pos, rho_b_ed_pos, rho_u_ed_neg, rho_b_ed_neg, rho_u_sc_pos, rho_b_sc_pos, rho_u_sc_neg, rho_b_sc_neg, rho_d_ed, rho_d_sc]
N_sl: [12]
b_sl: [2.86e-10]

11
examples/config/phase/mechanical/plastic/nonlocal_Ni.yaml
View File

@ -1,10 +1,13 @@
type: nonlocal
references:
C. Kords,
On the role of dislocation transport in the constitutive description of crystal plasticity,
RWTH Aachen 2013,
http://publications.rwth-aachen.de/record/229993/files/4862.pdf
- C. Kords,
On the role of dislocation transport in the constitutive description of crystal plasticity,
RWTH Aachen 2013,
http://publications.rwth-aachen.de/record/229993/files/4862.pdf
output: [rho_u_ed_pos, rho_b_ed_pos, rho_u_ed_neg, rho_b_ed_neg, rho_u_sc_pos, rho_b_sc_pos, rho_u_sc_neg, rho_b_sc_neg, rho_d_ed, rho_d_sc]
N_sl: [12]
b_sl: [2.48e-10]

4
examples/config/phase/mechanical/plastic/phenopowerlaw_Al.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- W.F. Hosford et al.,
Acta Metallurgica 8(3):187-199, 1960,
@ -7,8 +8,11 @@ references:
- U.F. Kocks,
Metallurgical and Materials Transactions B 1:1121–1143, 1970,
https://doi.org/10.1007/BF02900224
output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 20
a_sl: 3.1
h_0_sl-sl: 1.7e+8

4
examples/config/phase/mechanical/plastic/phenopowerlaw_Au.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- D. Ma et al.,
Acta Materialia 103:796-808, 2016,
@ -9,8 +10,11 @@ references:
- U.F. Kocks,
Metallurgical and Materials Transactions B 1:1121–1143, 1970,
https://doi.org/10.1007/BF02900224
output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 83.3
a_sl: 1.0
h_0_sl-sl: 75.0e+6

4
examples/config/phase/mechanical/plastic/phenopowerlaw_Cu.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- T Takeuchi,
Transactions of the Japan Institute of Metals 16(10):629-640, 1975,
@ -7,8 +8,11 @@ references:
- U.F. Kocks,
Metallurgical and Materials Transactions B 1:1121–1143, 1970,
https://doi.org/10.1007/BF02900224
output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 20
a_sl: 1.0
h_0_sl-sl: 2.4e+8

4
examples/config/phase/mechanical/plastic/phenopowerlaw_DP-steel-ferrite.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- C.C. Tasan et al.,
Acta Materialia 81:386-400, 2014,
@ -6,8 +7,11 @@ references:
- U.F. Kocks,
Metallurgical and Materials Transactions B 1:1121–1143, 1970,
https://doi.org/10.1007/BF02900224
output: [xi_sl, gamma_sl]
N_sl: [12, 12]
n_sl: 20
a_sl: 2.25
h_0_sl-sl: 1.0e+9

1
examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- F. Wang et al.,
Acta Materialia 80:77-93, 2014,

33
examples/config/phase/mechanical/plastic/phenopowerlaw_Sn-beta.yaml Normal file
View File

@ -0,0 +1,33 @@
type: phenopowerlaw
references:
- A. Chakraborty and P. Eisenlohr,
Journal of Applied Physics 124:025302, 2018,
https://doi.org/10.1063/1.5029933
output: [xi_sl, gamma_sl]
N_sl: [2, 2, 2, 4, 2, 4, 2, 2, 4, 0, 0, 8]
n_sl: 6.0
a_sl: 2.0
h_0_sl-sl: 20.0e+6
xi_0_sl: [8.5e+6, 4.3e+6, 10.4e+6, 4.5e+6, 5.6e+6, 5.1e+6, 7.4e+6, 15.0e+6, 6.6e+6, 0.0, 0.0, 12.0e+6]
xi_inf_sl: [11.0e+6, 9.0e+6, 11.0e+6, 9.0e+6, 10.0e+6, 10.0e+6, 10.0e+6, 10.0e+6, 9.0e+6, 0.0, 0.0, 13.0e+6]
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, # 50
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, # 100
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, # 150
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # unused entries are indicated by -1.0
dot_gamma_0_sl: 2.6e-8

3
examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- C. Zambaldi et al.,
Journal of Materials Research 27(1):356-367, 2021,
@ -6,9 +7,11 @@ references:
- L. Wang et al.,
Acta Materialia 132:598-610, 2017,
https://doi.org/10.1016/j.actamat.2017.05.015
output: [gamma_sl]
N_sl: [3, 3, 0, 12] # basal, prism, -, 1. pyr<c+a>
n_sl: 20
a_sl: 2.0
dot_gamma_0_sl: 0.001

1
examples/config/phase/thermal/AISI304.yaml
View File

@ -5,5 +5,6 @@ references:
- R.H. Bogaard et al.
Thermochimica Acta 218:373-393, 1993
https://doi.org/10.1016/0040-6031(93)80437-F
C_p: 470.0
K_11: 14.34

13
examples/config/phase/thermal/Al.yaml
View File

@ -1,5 +1,14 @@
references:
- https://www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
- J.G. Hust and A.B. Lankford,
Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
US Department of Commerce, Boulder, Colorado, 1984,
fit to Tab. 3.4.1 (RRR=1000, T_min=200K, T_max=900K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
K_11: 2.380e+2
K_11,T: 2.068e-3
K_11,T^2: -7.765e-5
T_ref: 293.15
C_p: 910.0
K_11: 236.0

1
examples/config/phase/thermal/Au.yaml
View File

@ -1,4 +1,5 @@
references:
- https://de.wikipedia.org/wiki/Gold
C_p: 128.0
K_11: 320.0

12
examples/config/phase/thermal/Cu.yaml
View File

@ -1,4 +1,14 @@
references:
- J.G. Hust and A.B. Lankford,
Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
US Department of Commerce, Boulder, Colorado, 1984,
fit to Tab. 2.4.1 (RRR=1000, T_min=200K, T_max=1000K)
- https://www.mit.edu/~6.777/matprops/copper.htm
K_11: 4.039e+2
K_11,T: -8.119e-2
K_11,T^2: 1.454e-5
T_ref: 293.15
C_p: 385.0
K_11: 401.0

14
examples/config/phase/thermal/Fe.yaml Normal file
View File

@ -0,0 +1,14 @@
references:
- J.G. Hust and A.B. Lankford,
Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
US Department of Commerce, Boulder, Colorado, 1984,
fit to Tab. 4.4.1 (RRR=300, T_min=200K, T_max=1000K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
K_11: 8.055e+1
K_11,T: -1.051e-1
K_11,T^2: 5.464e-5
T_ref: 293.15
C_p: 450.0

14
examples/config/phase/thermal/Ni.yaml Normal file
View File

@ -0,0 +1,14 @@
references:
- Y.S. Touloukian et al.,
TPRC Data Series Volume 1. Thermal conductivity - metallic elements and alloys,
IFI/Plenum, 1970,
fit to Tab. 35R (T_min=150K, T_max=500K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
K_11: 9.132e+1
K_11,T: -1.525e-1
K_11,T^2: 3.053e-4
T_ref: 293.15
C_p: 440.0

18
examples/config/phase/thermal/Sn-beta.yaml Normal file
View File

@ -0,0 +1,18 @@
references:
- Y.S. Touloukian et al.,
TPRC Data Series Volume 1. Thermal conductivity - metallic elements and alloys,
IFI/Plenum, 1970,
fit to Tab. 61R (T_min=100K, T_max=400K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
K_11: 7.414e+1
K_11,T: -6.465e-2
K_11,T^2: 2.066e-4
K_33: 5.147e+1
K_33,T: -4.506e-2
K_33,T^2: 1.435e-4
T_ref: 293.15
C_p: 210.0

12
examples/config/phase/thermal/W.yaml
View File

@ -1,4 +1,14 @@
references:
- J.G. Hust and A.B. Lankford,
Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
US Department of Commerce, Boulder, Colorado, 1984,
fit to Tab. 5.4.1 (RRR=300, T_min=200K, T_max=1000K)
- https://www.mit.edu/~6.777/matprops/tungsten.htm
K_11: 1.758e+2
K_11,T: -1.605e-1
K_11,T^2: 1.160e-4
T_ref: 293.15
C_p: 132.51
K_11: 178.0

1
examples/config/phase/thermal/source/dissipation_generic.yaml
View File

@ -1,2 +1,3 @@
type: dissipation
kappa: .9

1
examples/config/phase/thermal/source/externalheat_ramp-and-hold.yaml
View File

@ -1,3 +1,4 @@
type: externalheat
f_T: [1, 1, 0, 0]
t_n: [0, 500, 500.001, 1000]

1
examples/config/phase/thermal/steel-0.5C.yaml
View File

@ -1,5 +1,6 @@
references:
- https://www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
C_p: 490.0
K_11: 54.0

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha5-603-ge0ed668ce
v3.0.0-alpha5-696-g6fce27dee

22
python/damask/_configmaterial.py
View File

@ -178,8 +178,9 @@ class ConfigMaterial(Config):
table : damask.Table
Table that contains material information.
**kwargs
Keyword arguments where the key is the name and the value specifies
the label of the data column in the table.
Keyword arguments where the key is the property name and
the value specifies either the label of the data column in the table
or a constant value.
Returns
-------
@ -211,8 +212,23 @@ class ConfigMaterial(Config):
homogenization: {}
phase: {}
>>> cm.from_table(t,O='qu',phase='phase',homogenization='single_crystal')
material:
- constituents:
- O: [0.19, 0.8, 0.24, -0.51]
v: 1.0
phase: Aluminum
homogenization: single_crystal
- constituents:
- O: [0.8, 0.19, 0.24, -0.51]
v: 1.0
phase: Steel
homogenization: single_crystal
homogenization: {}
phase: {}
"""
kwargs_ = {k:table.get(v) for k,v in kwargs.items()}
kwargs_ = {k:table.get(v) if v in table.labels else np.atleast_2d([v]*len(table)).T for k,v in kwargs.items()}
_,idx = np.unique(np.hstack(list(kwargs_.values())),return_index=True,axis=0)
idx = np.sort(idx)

40
python/damask/_crystal.py
View File

@ -2,10 +2,11 @@ from typing import Union, Dict, List, Tuple
import numpy as np
from ._typehints import FloatSequence, CrystalFamily, CrystalLattice, CrystalKinematics
from . import util
from . import Rotation
lattice_symmetries = {
lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
'aP': 'triclinic',
'mP': 'monoclinic',
@ -30,9 +31,9 @@ lattice_symmetries = {
class Crystal():
"""Crystal lattice."""
def __init__(self,*,
family = None,
lattice = None,
def __init__(self, *,
family: CrystalFamily = None,
lattice: CrystalLattice = None,
a: float = None, b: float = None, c: float = None,
alpha: float = None, beta: float = None, gamma: float = None,
degrees: bool = False):
@ -130,9 +131,8 @@ class Crystal():
Crystal to check for equality.
"""
if not isinstance(other, Crystal):
return NotImplemented
return self.lattice == other.lattice and \
return NotImplemented if not isinstance(other, Crystal) else \
self.lattice == other.lattice and \
self.parameters == other.parameters and \
self.family == other.family
@ -208,7 +208,7 @@ class Crystal():
... }
"""
_basis = {
_basis: Dict[CrystalFamily, Dict[str, np.ndarray]] = {
'cubic': {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
@ -315,19 +315,19 @@ class Crystal():
self.lattice[-1],None),dtype=float)
def to_lattice(self, *,
direction: np.ndarray = None,
plane: np.ndarray = None) -> np.ndarray:
direction: FloatSequence = None,
plane: FloatSequence = None) -> np.ndarray:
"""
Calculate lattice vector corresponding to crystal frame direction or plane normal.
Parameters
----------
direction|plane : numpy.ndarray of shape (...,3)
direction|plane : numpy.ndarray, shape (...,3)
Vector along direction or plane normal.
Returns
-------
Miller : numpy.ndarray of shape (...,3)
Miller : numpy.ndarray, shape (...,3)
Lattice vector of direction or plane.
Use util.scale_to_coprime to convert to (integer) Miller indices.
@ -341,19 +341,19 @@ class Crystal():
def to_frame(self, *,
uvw: np.ndarray = None,
hkl: np.ndarray = None) -> np.ndarray:
uvw: FloatSequence = None,
hkl: FloatSequence = None) -> np.ndarray:
"""
Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl).
Parameters
----------
uvw|hkl : numpy.ndarray of shape (...,3)
uvw|hkl : numpy.ndarray, shape (...,3)
Miller indices of crystallographic direction or plane normal.
Returns
-------
vector : numpy.ndarray of shape (...,3)
vector : numpy.ndarray, shape (...,3)
Crystal frame vector along [uvw] direction or (hkl) plane normal.
"""
@ -366,7 +366,7 @@ class Crystal():
def kinematics(self,
mode: str) -> Dict[str, List[np.ndarray]]:
mode: CrystalKinematics) -> Dict[str, List[np.ndarray]]:
"""
Return crystal kinematics systems.
@ -381,7 +381,7 @@ class Crystal():
Directions and planes of deformation mode families.
"""
_kinematics = {
_kinematics: Dict[CrystalLattice, Dict[CrystalKinematics, List[np.ndarray]]] = {
'cF': {
'slip': [np.array([
[+0,+1,-1, +1,+1,+1],
@ -626,7 +626,7 @@ class Crystal():
def relation_operations(self,
model: str) -> Tuple[str, Rotation]:
model: str) -> Tuple[CrystalLattice, Rotation]:
"""
Crystallographic orientation relationships for phase transformations.
@ -658,7 +658,7 @@ class Crystal():
https://doi.org/10.1016/j.actamat.2004.11.021
"""
_orientation_relationships = {
_orientation_relationships: Dict[str, Dict[CrystalLattice,np.ndarray]] = {
'KS': {
'cF' : np.array([
[[-1, 0, 1],[ 1, 1, 1]],

238
python/damask/_orientation.py
View File

@ -1,8 +1,10 @@
import inspect
import copy
from typing import Union, Callable, List, Dict, Any, Tuple, TypeVar
import numpy as np
from ._typehints import FloatSequence, IntSequence, CrystalFamily, CrystalLattice
from . import Rotation
from . import Crystal
from . import util
@ -33,6 +35,7 @@ _parameter_doc = \
"""
MyType = TypeVar('MyType', bound='Orientation')
class Orientation(Rotation,Crystal):
"""
@ -93,12 +96,13 @@ class Orientation(Rotation,Crystal):
@util.extend_docstring(_parameter_doc)
def __init__(self,
rotation = np.array([1.0,0.0,0.0,0.0]), *,
family = None,
lattice = None,
a = None,b = None,c = None,
alpha = None,beta = None,gamma = None,
degrees = False):
rotation: Union[FloatSequence, Rotation] = np.array([1.,0.,0.,0.]),
*,
family: CrystalFamily = None,
lattice: CrystalLattice = None,
a: float = None, b: float = None, c: float = None,
alpha: float = None, beta: float = None, gamma: float = None,
degrees: bool = False):
"""
New orientation.
@ -115,13 +119,13 @@ class Orientation(Rotation,Crystal):
a=a,b=b,c=c, alpha=alpha,beta=beta,gamma=gamma, degrees=degrees)
def __repr__(self):
def __repr__(self) -> str:
"""Represent."""
return '\n'.join([Crystal.__repr__(self),
Rotation.__repr__(self)])
def __copy__(self,rotation=None):
def __copy__(self: MyType,
rotation: Union[FloatSequence, Rotation] = None) -> MyType:
"""Create deep copy."""
dup = copy.deepcopy(self)
if rotation is not None:
@ -131,7 +135,9 @@ class Orientation(Rotation,Crystal):
copy = __copy__
def __eq__(self,other):
def __eq__(self,
other: object) -> bool:
"""
Equal to other.
@ -141,12 +147,15 @@ class Orientation(Rotation,Crystal):
Orientation to check for equality.
"""
if not isinstance(other, Orientation):
return NotImplemented
matching_type = self.family == other.family and \
self.lattice == other.lattice and \
self.parameters == other.parameters
return np.logical_and(matching_type,super(self.__class__,self.reduced).__eq__(other.reduced))
def __ne__(self,other):
def __ne__(self,
other: object) -> bool:
"""
Not equal to other.
@ -156,10 +165,14 @@ class Orientation(Rotation,Crystal):
Orientation to check for equality.
"""
return np.logical_not(self==other)
return np.logical_not(self==other) if isinstance(other, Orientation) else NotImplemented
def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
def isclose(self: MyType,
other: MyType,
rtol: float = 1e-5,
atol: float = 1e-8,
equal_nan: bool = True) -> bool:
"""
Report where values are approximately equal to corresponding ones of other Orientation.
@ -176,7 +189,7 @@ class Orientation(Rotation,Crystal):
Returns
-------
mask : numpy.ndarray bool
mask : numpy.ndarray of bool, shape (self.shape)
Mask indicating where corresponding orientations are close.
"""
@ -187,7 +200,11 @@ class Orientation(Rotation,Crystal):
def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
def allclose(self: MyType,
other: MyType,
rtol: float = 1e-5,
atol: float = 1e-8,
equal_nan: bool = True) -> bool:
"""
Test whether all values are approximately equal to corresponding ones of other Orientation.
@ -208,10 +225,11 @@ class Orientation(Rotation,Crystal):
Whether all values are close between both orientations.
"""
return np.all(self.isclose(other,rtol,atol,equal_nan))
return bool(np.all(self.isclose(other,rtol,atol,equal_nan)))
def __mul__(self,other):
def __mul__(self: MyType,
other: Union[Rotation, 'Orientation']) -> MyType:
"""
Compose this orientation with other.
@ -226,14 +244,15 @@ class Orientation(Rotation,Crystal):
Compound rotation self*other, i.e. first other then self rotation.
"""
if isinstance(other,Orientation) or isinstance(other,Rotation):
return self.copy(rotation=Rotation.__mul__(self,Rotation(other.quaternion)))
if isinstance(other, (Orientation,Rotation)):
return self.copy(Rotation(self.quaternion)*Rotation(other.quaternion))
else:
raise TypeError('use "O@b", i.e. matmul, to apply Orientation "O" to object "b"')
@staticmethod
def _split_kwargs(kwargs,target):
def _split_kwargs(kwargs: Dict[str, Any],
target: Callable) -> Tuple[Dict[str, Any], ...]:
"""
Separate keyword arguments in 'kwargs' targeted at 'target' from general keyword arguments of Orientation objects.
@ -252,7 +271,7 @@ class Orientation(Rotation,Crystal):
Valid keyword arguments of Orientation object.
"""
kws = ()
kws: Tuple[Dict[str, Any], ...] = ()
for t in (target,Orientation.__init__):
kws += ({key: kwargs[key] for key in set(inspect.signature(t).parameters) & set(kwargs)},)
@ -264,105 +283,108 @@ class Orientation(Rotation,Crystal):
@classmethod
@util.extended_docstring(Rotation.from_random,_parameter_doc)
def from_random(cls,**kwargs):
@util.extended_docstring(Rotation.from_random, _parameter_doc)
def from_random(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_random)
return cls(rotation=Rotation.from_random(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_quaternion,_parameter_doc)
def from_quaternion(cls,**kwargs):
def from_quaternion(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_quaternion)
return cls(rotation=Rotation.from_quaternion(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc)
def from_Euler_angles(cls,**kwargs):
def from_Euler_angles(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Euler_angles)
return cls(rotation=Rotation.from_Euler_angles(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_axis_angle,_parameter_doc)
def from_axis_angle(cls,**kwargs):
def from_axis_angle(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_axis_angle)
return cls(rotation=Rotation.from_axis_angle(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_basis,_parameter_doc)
def from_basis(cls,**kwargs):
def from_basis(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_basis)
return cls(rotation=Rotation.from_basis(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_matrix,_parameter_doc)
def from_matrix(cls,**kwargs):
def from_matrix(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_matrix)
return cls(rotation=Rotation.from_matrix(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc)
def from_Rodrigues_vector(cls,**kwargs):
def from_Rodrigues_vector(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Rodrigues_vector)
return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_homochoric,_parameter_doc)
def from_homochoric(cls,**kwargs):
def from_homochoric(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_homochoric)
return cls(rotation=Rotation.from_homochoric(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_cubochoric,_parameter_doc)
def from_cubochoric(cls,**kwargs):
def from_cubochoric(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_cubochoric)
return cls(rotation=Rotation.from_cubochoric(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_spherical_component,_parameter_doc)
def from_spherical_component(cls,**kwargs):
def from_spherical_component(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_spherical_component)
return cls(rotation=Rotation.from_spherical_component(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_fiber_component,_parameter_doc)
def from_fiber_component(cls,**kwargs):
def from_fiber_component(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_fiber_component)
return cls(rotation=Rotation.from_fiber_component(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extend_docstring(_parameter_doc)
def from_directions(cls,uvw,hkl,**kwargs):
def from_directions(cls,
uvw: FloatSequence,
hkl: FloatSequence,
**kwargs) -> 'Orientation':
"""
Initialize orientation object from two crystallographic directions.
Parameters
----------
uvw : list, numpy.ndarray of shape (...,3)
lattice direction aligned with lab frame x-direction.
hkl : list, numpy.ndarray of shape (...,3)
lattice plane normal aligned with lab frame z-direction.
uvw : numpy.ndarray, shape (...,3)
Lattice direction aligned with lab frame x-direction.
hkl : numpy.ndarray, shape (...,3)
Lattice plane normal aligned with lab frame z-direction.
"""
o = cls(**kwargs)
x = o.to_frame(uvw=uvw)
z = o.to_frame(hkl=hkl)
om = np.stack([x,np.cross(z,x),z],axis=-2)
return o.copy(rotation=Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
return o.copy(Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
@property
def equivalent(self):
def equivalent(self: MyType) -> MyType:
"""
Orientations that are symmetrically equivalent.
@ -372,11 +394,11 @@ class Orientation(Rotation,Crystal):
"""
sym_ops = self.symmetry_operations
o = sym_ops.broadcast_to(sym_ops.shape+self.shape,mode='right')
return self.copy(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
return self.copy(o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
@property
def reduced(self):
def reduced(self: MyType) -> MyType:
"""Select symmetrically equivalent orientation that falls into fundamental zone according to symmetry."""
eq = self.equivalent
ok = eq.in_FZ
@ -387,13 +409,13 @@ class Orientation(Rotation,Crystal):
@property
def in_FZ(self):
def in_FZ(self) -> Union[np.bool_, np.ndarray]:
"""
Check whether orientation falls into fundamental zone of own symmetry.
Returns
-------
in : numpy.ndarray of bool, quaternion.shape
in : numpy.ndarray of bool, shape (self.shape)
Whether Rodrigues-Frank vector falls into fundamental zone.
Notes
@ -431,13 +453,13 @@ class Orientation(Rotation,Crystal):
@property
def in_disorientation_FZ(self):
def in_disorientation_FZ(self) -> np.ndarray:
"""
Check whether orientation falls into fundamental zone of disorientations.
Returns
-------
in : numpy.ndarray of bool, quaternion.shape
in : numpy.ndarray of bool, shape (self.shape)
Whether Rodrigues-Frank vector falls into disorientation FZ.
References
@ -471,8 +493,9 @@ class Orientation(Rotation,Crystal):
else:
return np.ones_like(rho[...,0],dtype=bool)
def disorientation(self,other,return_operators=False):
def disorientation(self,
other: 'Orientation',
return_operators: bool = False) -> object:
"""
Calculate disorientation between myself and given other orientation.
@ -490,7 +513,7 @@ class Orientation(Rotation,Crystal):
-------
disorientation : Orientation
Disorientation between self and other.
operators : numpy.ndarray int of shape (...,2), conditional
operators : numpy.ndarray of int, shape (...,2), conditional
Index of symmetrically equivalent orientation that rotated vector to the SST.
Notes
@ -530,7 +553,7 @@ class Orientation(Rotation,Crystal):
raise NotImplementedError('disorientation between different crystal families')
blend = util.shapeblender(self.shape,other.shape)
s = self.equivalent
s = self.equivalent
o = other.equivalent
s_ = s.reshape((s.shape[0],1)+ self.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right')
@ -557,13 +580,15 @@ class Orientation(Rotation,Crystal):
)
def average(self,weights=None,return_cloud=False):
def average(self,
weights: FloatSequence = None,
return_cloud: bool = False):
"""
Return orientation average over last dimension.
Parameters
----------
weights : numpy.ndarray, optional
weights : numpy.ndarray, shape (self.shape), optional
Relative weights of orientations.
return_cloud : bool, optional
Return the set of symmetrically equivalent orientations that was used in averaging.
@ -583,31 +608,30 @@ class Orientation(Rotation,Crystal):
"""
eq = self.equivalent
m = eq.misorientation(self[...,0].reshape((1,)+self.shape[:-1]+(1,))
.broadcast_to(eq.shape))\
.as_axis_angle()[...,3]
m = eq.misorientation(self[...,0].reshape((1,)+self.shape[:-1]+(1,)) # type: ignore
.broadcast_to(eq.shape)).as_axis_angle()[...,3] # type: ignore
r = Rotation(np.squeeze(np.take_along_axis(eq.quaternion,
np.argmin(m,axis=0)[np.newaxis,...,np.newaxis],
axis=0),
axis=0))
return (
(self.copy(rotation=Rotation(r).average(weights)),
self.copy(rotation=Rotation(r)))
if return_cloud else
self.copy(rotation=Rotation(r).average(weights))
return ((self.copy(Rotation(r).average(weights)),self.copy(Rotation(r))) if return_cloud else
self.copy(Rotation(r).average(weights))
)
def to_SST(self,vector,proper=False,return_operators=False):
def to_SST(self,
vector: FloatSequence,
proper: bool = False,
return_operators: bool = False) -> np.ndarray:
"""
Rotate vector to ensure it falls into (improper or proper) standard stereographic triangle of crystal symmetry.
Parameters
----------
vector : numpy.ndarray of shape (...,3)
vector : numpy.ndarray, shape (...,3)
Lab frame vector to align with crystal frame direction.
Shape of vector blends with shape of own rotation array.
For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs.
For example, a rotation array of shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
proper : bool, optional
Consider only vectors with z >= 0, hence combine two neighboring SSTs.
Defaults to False.
@ -617,15 +641,18 @@ class Orientation(Rotation,Crystal):
Returns
-------
vector_SST : numpy.ndarray of shape (...,3)
vector_SST : numpy.ndarray, shape (...,3)
Rotated vector falling into SST.
operators : numpy.ndarray int of shape (...), conditional
operators : numpy.ndarray of int, shape (...), conditional
Index of symmetrically equivalent orientation that rotated vector to SST.
"""
vector_ = np.array(vector,float)
if vector_.shape[-1] != 3:
raise ValueError('input is not a field of three-dimensional vectors')
eq = self.equivalent
blend = util.shapeblender(eq.shape,np.array(vector).shape[:-1])
poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(np.array(vector),blend+(3,))
blend = util.shapeblender(eq.shape,vector_.shape[:-1])
poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(vector_,blend+(3,))
ok = self.in_SST(poles,proper=proper)
ok &= np.cumsum(ok,axis=0) == 1
loc = np.where(ok)
@ -637,13 +664,15 @@ class Orientation(Rotation,Crystal):
)
def in_SST(self,vector,proper=False):
def in_SST(self,
vector: FloatSequence,
proper: bool = False) -> Union[np.bool_, np.ndarray]:
"""
Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry.
Parameters
----------
vector : numpy.ndarray of shape (...,3)
vector : numpy.ndarray, shape (...,3)
Vector to check.
proper : bool, optional
Consider only vectors with z >= 0, hence combine two neighboring SSTs.
@ -655,39 +684,43 @@ class Orientation(Rotation,Crystal):
Whether vector falls into SST.
"""
if not isinstance(vector,np.ndarray) or vector.shape[-1] != 3:
vector_ = np.array(vector,float)
if vector_.shape[-1] != 3:
raise ValueError('input is not a field of three-dimensional vectors')
if self.standard_triangle is None: # direct exit for no symmetry
return np.ones_like(vector[...,0],bool)
return np.ones_like(vector_[...,0],bool)
if proper:
components_proper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['proper'], vector.shape+(3,)),
vector), 12)
np.broadcast_to(self.standard_triangle['proper'], vector_.shape+(3,)),
vector_), 12)
components_improper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)),
vector), 12)
np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
vector_), 12)
return np.all(components_proper >= 0.0,axis=-1) \
| np.all(components_improper >= 0.0,axis=-1)
else:
components = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)),
np.block([vector[...,:2],np.abs(vector[...,2:3])])), 12)
np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
return np.all(components >= 0.0,axis=-1)
def IPF_color(self,vector,in_SST=True,proper=False):
def IPF_color(self,
vector: FloatSequence,
in_SST: bool = True,
proper: bool = False) -> np.ndarray:
"""
Map vector to RGB color within standard stereographic triangle of own symmetry.
Parameters
----------
vector : numpy.ndarray of shape (...,3)
vector : numpy.ndarray, shape (...,3)
Vector to colorize.
Shape of vector blends with shape of own rotation array.
For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs.
For example, a rotation array of shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
in_SST : bool, optional
Consider symmetrically equivalent orientations such that poles are located in SST.
Defaults to True.
@ -697,7 +730,7 @@ class Orientation(Rotation,Crystal):
Returns
-------
rgb : numpy.ndarray of shape (...,3)
rgb : numpy.ndarray, shape (...,3)
RGB array of IPF colors.
Examples
@ -721,35 +754,35 @@ class Orientation(Rotation,Crystal):
if proper:
components_proper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['proper'], vector_.shape+(3,)),
np.broadcast_to(self.standard_triangle['proper'], vector_.shape+(3,)),
vector_), 12)
components_improper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
vector_), 12)
in_SST = np.all(components_proper >= 0.0,axis=-1) \
| np.all(components_improper >= 0.0,axis=-1)
components = np.where((in_SST & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis],
in_SST_ = np.all(components_proper >= 0.0,axis=-1) \
| np.all(components_improper >= 0.0,axis=-1)
components = np.where((in_SST_ & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis],
components_proper,components_improper)
else:
components = np.around(np.einsum('...ji,...i',
np.broadcast_to(self .standard_triangle['improper'], vector_.shape+(3,)),
np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
in_SST = np.all(components >= 0.0,axis=-1)
in_SST_ = np.all(components >= 0.0,axis=-1)
with np.errstate(invalid='ignore',divide='ignore'):
rgb = (components/np.linalg.norm(components,axis=-1,keepdims=True))**0.5 # smoothen color ramps
rgb = np.clip(rgb,0.,1.) # clip intensity
rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1
rgb[np.broadcast_to(~in_SST[...,np.newaxis],rgb.shape)] = 0.0
rgb[np.broadcast_to(~in_SST_[...,np.newaxis],rgb.shape)] = 0.0
return rgb
@property
def symmetry_operations(self):
def symmetry_operations(self) -> Rotation:
"""Symmetry operations as Rotations."""
_symmetry_operations = {
_symmetry_operations: Dict[CrystalFamily, List] = {
'cubic': [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
@ -819,22 +852,25 @@ class Orientation(Rotation,Crystal):
####################################################################################################
# functions that require lattice, not just family
def to_pole(self,*,uvw=None,hkl=None,with_symmetry=False):
def to_pole(self, *,
uvw: FloatSequence = None,
hkl: FloatSequence = None,
with_symmetry: bool = False) -> np.ndarray:
"""
Calculate lab frame vector along lattice direction [uvw] or plane normal (hkl).
Parameters
----------
uvw|hkl : numpy.ndarray of shape (...,3)
uvw|hkl : numpy.ndarray, shape (...,3)
Miller indices of crystallographic direction or plane normal.
Shape of vector blends with shape of own rotation array.
For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs.
For example, a rotation array, shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
with_symmetry : bool, optional
Calculate all N symmetrically equivalent vectors.
Returns
-------
vector : numpy.ndarray of shape (...,3) or (...,N,3)
vector : numpy.ndarray, shape (...,3) or (...,N,3)
Lab frame vector (or vectors if with_symmetry) along [uvw] direction or (hkl) plane normal.
"""
@ -846,23 +882,24 @@ class Orientation(Rotation,Crystal):
blend += sym_ops.shape
v = sym_ops.broadcast_to(shape) \
@ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,))
return ~(self.broadcast_to(blend)) \
@ np.broadcast_to(v,blend+(3,))
return ~(self.broadcast_to(blend))@ np.broadcast_to(v,blend+(3,))
def Schmid(self,*,N_slip=None,N_twin=None):
def Schmid(self, *,
N_slip: IntSequence = None,
N_twin: IntSequence = None) -> np.ndarray:
u"""
Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems.
Parameters
----------
N_slip|N_twin : iterable of int
N_slip|N_twin : '*' or iterable of int
Number of deformation systems per family of the deformation system.
Use '*' to select all.
Returns
-------
P : numpy.ndarray of shape (N,...,3,3)
P : numpy.ndarray, shape (N,...,3,3)
Schmid matrix for each of the N deformation systems.
Examples
@ -887,6 +924,8 @@ class Orientation(Rotation,Crystal):
(self.kinematics('twin'),N_twin)
if active == '*': active = [len(a) for a in kinematics['direction']]
if not active:
raise ValueError('Schmid matrix not defined')
d = self.to_frame(uvw=np.vstack([kinematics['direction'][i][:n] for i,n in enumerate(active)]))
p = self.to_frame(hkl=np.vstack([kinematics['plane'][i][:n] for i,n in enumerate(active)]))
P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=1,keepdims=True),
@ -897,7 +936,8 @@ class Orientation(Rotation,Crystal):
@ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape)
def related(self,model):
def related(self: MyType,
model: str) -> MyType:
"""
Orientations derived from the given relationship.

418
python/damask/_rotation.py
View File

File diff suppressed because it is too large Load Diff

39
python/damask/_table.py
View File

@ -44,6 +44,13 @@ class Table:
return '\n'.join(['# '+c for c in self.comments])+'\n'+data_repr
def __eq__(self,
other: object) -> bool:
"""Compare to other Table."""
return NotImplemented if not isinstance(other,Table) else \
self.shapes == other.shapes and self.data.equals(other.data)
def __getitem__(self,
item: Union[slice, Tuple[slice, ...]]) -> 'Table':
"""
@ -75,20 +82,22 @@ class Table:
colB colA
0 1 0
2 7 6
>>> tbl[1:2,'colB']
>>> tbl[[True,False,False,True],'colB']
colB
1 4
2 7
0 1
3 10
"""
item = (item,slice(None,None,None)) if isinstance(item,slice) else \
item if isinstance(item[0],slice) else \
(slice(None,None,None),item)
sliced = self.data.loc[item]
cols = np.array(sliced.columns if isinstance(sliced,pd.core.frame.DataFrame) else [item[1]])
item_ = (item,slice(None,None,None)) if isinstance(item,(slice,np.ndarray)) else \
(np.array(item),slice(None,None,None)) if isinstance(item,list) and np.array(item).dtype == np.bool_ else \
(np.array(item[0]),item[1]) if isinstance(item[0],list) else \
item if isinstance(item[0],(slice,np.ndarray)) else \
(slice(None,None,None),item)
sliced = self.data.loc[item_]
cols = np.array(sliced.columns if isinstance(sliced,pd.core.frame.DataFrame) else [item_[1]])
_,idx = np.unique(cols,return_index=True)
return self.__class__(data=sliced,
shapes = {k:self.shapes[k] for k in cols[np.sort(idx)]},
shapes={k:self.shapes[k] for k in cols[np.sort(idx)]},
comments=self.comments)
@ -185,7 +194,7 @@ class Table:
Returns
-------
mask : numpy.ndarray bool
mask : numpy.ndarray of bool
Mask indicating where corresponding table values are close.
"""
@ -361,9 +370,9 @@ class Table:
label : str
Column label.
data : numpy.ndarray
New data.
Replacement data.
info : str, optional
Human-readable information about the new data.
Human-readable information about the modified data.
Returns
-------
@ -396,9 +405,9 @@ class Table:
label : str
Column label.
data : numpy.ndarray
Modified data.
New data.
info : str, optional
Human-readable information about the modified data.
Human-readable information about the new data.
Returns
-------
@ -474,7 +483,7 @@ class Table:
labels: Union[str, List[str]],
ascending: Union[bool, List[bool]] = True) -> 'Table':
"""
Sort table by values of given labels.
Sort table by data of given columns.
Parameters
----------

5
python/damask/_typehints.py
View File

@ -1,6 +1,6 @@
"""Functionality for typehints."""
from typing import Sequence, Union, TextIO
from typing import Sequence, Union, Literal, TextIO
from pathlib import Path
import numpy as np
@ -9,6 +9,9 @@ import numpy as np
FloatSequence = Union[np.ndarray,Sequence[float]]
IntSequence = Union[np.ndarray,Sequence[int]]
FileHandle = Union[TextIO, str, Path]
CrystalFamily = Union[None,Literal['triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic']]
CrystalLattice = Union[None,Literal['aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF']]
CrystalKinematics = Literal['slip', 'twin']
NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.Generator]
# BitGenerator does not exists in older numpy versions
#NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.BitGenerator, np.random.Generator]

6
python/damask/util.py
View File

@ -9,7 +9,7 @@ import re
import fractions
from collections import abc
from functools import reduce
from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal, SupportsIndex, Sequence
from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal
from pathlib import Path
import numpy as np
@ -427,7 +427,7 @@ def hybrid_IA(dist: np.ndarray,
def shapeshifter(fro: Tuple[int, ...],
to: Tuple[int, ...],
mode: Literal['left','right'] = 'left',
keep_ones: bool = False) -> Sequence[SupportsIndex]:
keep_ones: bool = False) -> Tuple[int, ...]:
"""
Return dimensions that reshape 'fro' to become broadcastable to 'to'.
@ -490,7 +490,7 @@ def shapeshifter(fro: Tuple[int, ...],
def shapeblender(a: Tuple[int, ...],
b: Tuple[int, ...]) -> Sequence[SupportsIndex]:
b: Tuple[int, ...]) -> Tuple[int, ...]:
"""
Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.

12
python/tests/test_ConfigMaterial.py
View File

@ -96,6 +96,18 @@ class TestConfigMaterial:
for i,m in enumerate(c['material']):
assert m['homogenization'] == 1 and (m['constituents'][0]['O'] == [1,0,1,1]).all()
def test_from_table_with_constant(self):
N = np.random.randint(3,10)
a = np.vstack((np.hstack((np.arange(N),np.arange(N)[::-1])),
np.ones(N*2),np.zeros(N*2),np.ones(N*2),np.ones(N*2),
np.ones(N*2),
)).T
t = Table(a,{'varying':1,'constant':4,'ones':1})
c = ConfigMaterial.from_table(t,**{'phase':'varying','O':'constant','homogenization':1})
assert len(c['material']) == N
for i,m in enumerate(c['material']):
assert m['homogenization'] == 1 and (m['constituents'][0]['O'] == [1,0,1,1]).all()
@pytest.mark.parametrize('N,n,kw',[
(1,1,{'phase':'Gold',
'O':[1,0,0,0],

6
python/tests/test_Table.py
View File

@ -59,10 +59,14 @@ class TestTable:
@pytest.mark.parametrize('N',[1,3,4])
def test_slice(self,default,N):
mask = np.random.choice([True,False],len(default))
assert len(default[:N]) == 1+N
assert len(default[:N,['F','s']]) == 1+N
assert len(default[mask,['F','s']]) == np.count_nonzero(mask)
assert default[mask,['F','s']] == default[mask][['F','s']] == default[['F','s']][mask]
assert default[np.logical_not(mask),['F','s']] != default[mask][['F','s']]
assert default[N:].get('F').shape == (len(default)-N,3,3)
assert (default[:N,['v','s']].data == default['v','s'][:N].data).all().all()
assert default[:N,['v','s']].data.equals(default['v','s'][:N].data)
@pytest.mark.parametrize('mode',['str','path'])
def test_write_read(self,default,tmp_path,mode):

3
src/CPFEM.f90
View File

@ -192,11 +192,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else validCalculation
if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP])
call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip)
if (.not. terminallyIll) &
call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP])
terminalIllness: if (terminallyIll) then
call random_number(rnd)

4
src/discretization.f90
View File

@ -12,7 +12,8 @@ module discretization
integer, public, protected :: &
discretization_nIPs, &
discretization_Nelems
discretization_Nelems, &
discretization_Ncells
integer, public, protected, dimension(:), allocatable :: &
discretization_materialAt !ToDo: discretization_ID_material
@ -53,6 +54,7 @@ subroutine discretization_init(materialAt,&
discretization_Nelems = size(materialAt,1)
discretization_nIPs = size(IPcoords0,2)/discretization_Nelems
discretization_Ncells = discretization_Nelems*discretization_nIPs
discretization_materialAt = materialAt

182
src/grid/spectral_utilities.f90
View File

@ -31,7 +31,7 @@ module spectral_utilities
!--------------------------------------------------------------------------------------------------
! grid related information
real(pReal), protected, public :: wgt !< weighting factor 1/Nelems
integer, protected, public :: grid1Red !< cells(1)/2
integer, protected, public :: cells1Red !< cells(1)/2
real(pReal), protected, public, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence
!--------------------------------------------------------------------------------------------------
@ -201,7 +201,7 @@ subroutine spectral_utilities_init
num_grid%get_asString('PETSc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc)
grid1Red = cells(1)/2 + 1
cells1Red = cells(1)/2 + 1
wgt = 1.0/real(product(cells),pReal)
num%memory_efficient = num_grid%get_asInt('memory_efficient', defaultVal=1) > 0 ! ToDo: should be logical in YAML file
@ -265,8 +265,8 @@ subroutine spectral_utilities_init
gridFFTW = int(cells,C_INTPTR_T)
alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, &
PETSC_COMM_WORLD, local_K, local_K_offset)
allocate (xi1st (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension
allocate (xi2nd (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension
allocate (xi1st (3,cells1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension
allocate (xi2nd (3,cells1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension
tensorField = fftw_alloc_complex(tensorSize*alloc_local)
call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, &
@ -333,7 +333,7 @@ subroutine spectral_utilities_init
do j = 1, cells(2)
k_s(2) = j - 1
if (j > cells(2)/2 + 1) k_s(2) = k_s(2) - cells(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do i = 1, grid1Red
do i = 1, cells1Red
k_s(1) = i - 1 ! symmetry, junst running from 0,1,...,N/2,N/2+1
xi2nd(1:3,i,j,k-cells3Offset) = utilities_getFreqDerivative(k_s)
where(mod(cells,2)==0 .and. [i,j,k] == cells/2+1 .and. &
@ -347,7 +347,7 @@ subroutine spectral_utilities_init
if (num%memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal))
else ! precalculation of gamma_hat field
allocate (gamma_hat(3,3,3,3,grid1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
allocate (gamma_hat(3,3,3,3,cells1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
endif
end subroutine spectral_utilities_init
@ -362,7 +362,7 @@ end subroutine spectral_utilities_init
subroutine utilities_updateGamma(C)
real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx
real(pReal), dimension(6,6) :: A, A_inv
integer :: &
i, j, k, &
@ -372,27 +372,40 @@ subroutine utilities_updateGamma(C)
C_ref = C
if (.not. num%memory_efficient) then
gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, grid1Red
gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A
!$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err)
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells1Red
if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
do concurrent (l = 1:3, m = 1:3)
#ifndef __INTEL_COMPILER
do concurrent(l = 1:3, m = 1:3)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
end do
do concurrent(l = 1:3, m = 1:3)
temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
end do
A(1:3,1:3) = temp33_complex%re; A(4:6,4:6) = temp33_complex%re
A(1:3,4:6) = temp33_complex%im; A(4:6,1:3) = -temp33_complex%im
#else
forall(l = 1:3, m = 1:3) &
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
forall(l = 1:3, m = 1:3) &
temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
#endif
A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) = temp33_cmplx%re
A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
call math_invert(A_inv, err, A)
temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
#ifndef __INTEL_COMPILER
do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_complex(l,n)* &
conjg(-xi1st(o,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
end do
#else
forall(l=1:3, m=1:3, n=1:3, o=1:3) &
gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
#endif
end if
end if
end do; end do; end do
!$OMP END PARALLEL DO
endif
end subroutine utilities_updateGamma
@ -405,7 +418,7 @@ end subroutine utilities_updateGamma
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorForward
tensorField_real(1:3,1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
end subroutine utilities_FFTtensorForward
@ -429,7 +442,7 @@ end subroutine utilities_FFTtensorBackward
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarForward
scalarField_real(cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
end subroutine utilities_FFTscalarForward
@ -454,7 +467,7 @@ end subroutine utilities_FFTscalarBackward
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorForward
vectorField_real(1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
end subroutine utilities_FFTvectorForward
@ -478,7 +491,7 @@ end subroutine utilities_FFTvectorBackward
subroutine utilities_fourierGammaConvolution(fieldAim)
real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx
real(pReal), dimension(6,6) :: A, A_inv
integer :: &
@ -493,38 +506,61 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium)
memoryEfficient: if (num%memory_efficient) then
do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red
!$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat)
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells1Red
if (any([i,j,k+cells3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
#ifndef __INTEL_COMPILER
do concurrent(l = 1:3, m = 1:3)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
end do
do concurrent(l = 1:3, m = 1:3)
temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
end do
A(1:3,1:3) = temp33_complex%re; A(4:6,4:6) = temp33_complex%re
A(1:3,4:6) = temp33_complex%im; A(4:6,1:3) = -temp33_complex%im
#else
forall(l = 1:3, m = 1:3) &
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
forall(l = 1:3, m = 1:3) &
temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
#endif
A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) = temp33_cmplx%re
A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
call math_invert(A_inv, err, A)
temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
#ifndef __INTEL_COMPILER
do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
gamma_hat(l,m,n,o,1,1,1) = temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k)
gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m)
end do
do concurrent(l = 1:3, m = 1:3)
temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
end do
#else
forall(l=1:3, m=1:3, n=1:3, o=1:3) &
gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m)
forall(l = 1:3, m = 1:3) &
temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
#endif
tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx
else
gamma_hat(1:3,1:3,1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal)
tensorField_fourier(1:3,1:3,i,j,k) = cmplx(0.0_pReal,0.0_pReal,pReal)
end if
do concurrent(l = 1:3, m = 1:3)
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
end do
tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
end if
end do; end do; end do
!$OMP END PARALLEL DO
else memoryEfficient
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
!$OMP PARALLEL DO PRIVATE(l,m,temp33_cmplx)
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells1Red
#ifndef __INTEL_COMPILER
do concurrent(l = 1:3, m = 1:3)
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k))
temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k))
end do
tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
#else
forall(l = 1:3, m = 1:3) &
temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k))
#endif
tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx
end do; end do; end do
!$OMP END PARALLEL DO
end if memoryEfficient
if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
@ -544,12 +580,14 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation
do k = 1, cells3; do j = 1, cells(2) ;do i = 1, grid1Red
!$OMP PARALLEL DO PRIVATE(GreenOp_hat)
do k = 1, cells3; do j = 1, cells(2) ;do i = 1, cells1Red
GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) &
/ (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal) &
* sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k))))
scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat
enddo; enddo; enddo
!$OMP END PARALLEL DO
end subroutine utilities_fourierGreenConvolution
@ -572,7 +610,7 @@ real(pReal) function utilities_divergenceRMS()
! calculating RMS divergence criterion in Fourier space
utilities_divergenceRMS = 0.0_pReal
do k = 1, cells3; do j = 1, cells(2)
do i = 2, grid1Red -1 ! Has somewhere a conj. complex counterpart. Therefore count it twice.
do i = 2, cells1Red -1 ! Has somewhere a conj. complex counterpart. Therefore count it twice.
utilities_divergenceRMS = utilities_divergenceRMS &
+ 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k), & ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2, i.e. do not take square root and square again
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2) & ! --> sum squared L_2 norm of vector
@ -584,10 +622,10 @@ real(pReal) function utilities_divergenceRMS()
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
+ sum( real(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2) &
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2)
+ sum( real(matmul(tensorField_fourier(1:3,1:3,cells1Red,j,k), &
conjg(-xi1st(1:3,cells1Red,j,k))*rescaledGeom))**2) &
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,cells1Red,j,k), &
conjg(-xi1st(1:3,cells1Red,j,k))*rescaledGeom))**2)
enddo; enddo
if (cells(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of cells(1) == 1
call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@ -617,7 +655,7 @@ real(pReal) function utilities_curlRMS()
utilities_curlRMS = 0.0_pReal
do k = 1, cells3; do j = 1, cells(2);
do i = 2, grid1Red - 1
do i = 2, cells1Red - 1
do l = 1, 3
curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3))
@ -640,12 +678,12 @@ real(pReal) function utilities_curlRMS()
utilities_curlRMS = utilities_curlRMS &
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (DC) does not have a conjugate complex counterpart (if cells(1) /= 1)
do l = 1, 3
curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3))
curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3) &
-tensorField_fourier(l,3,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1))
curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1) &
-tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2))
curl_fourier = (+tensorField_fourier(l,3,cells1Red,j,k)*xi1st(2,cells1Red,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,cells1Red,j,k)*xi1st(3,cells1Red,j,k)*rescaledGeom(3))
curl_fourier = (+tensorField_fourier(l,1,cells1Red,j,k)*xi1st(3,cells1Red,j,k)*rescaledGeom(3) &
-tensorField_fourier(l,3,cells1Red,j,k)*xi1st(1,cells1Red,j,k)*rescaledGeom(1))
curl_fourier = (+tensorField_fourier(l,2,cells1Red,j,k)*xi1st(1,cells1Red,j,k)*rescaledGeom(1) &
-tensorField_fourier(l,1,cells1Red,j,k)*xi1st(2,cells1Red,j,k)*rescaledGeom(2))
enddo
utilities_curlRMS = utilities_curlRMS &
+ sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (Nyquist) does not have a conjugate complex counterpart (if cells(1) /= 1)
@ -736,9 +774,10 @@ subroutine utilities_fourierScalarGradient()
integer :: i, j, k
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells1Red
vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) ! ToDo: no -conjg?
enddo; enddo; enddo
end do; end do; end do
end subroutine utilities_fourierScalarGradient
@ -748,11 +787,9 @@ end subroutine utilities_fourierScalarGradient
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierVectorDivergence()
integer :: i, j, k
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
scalarField_fourier(i,j,k) = sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k)))
enddo; enddo; enddo
scalarField_fourier(1:cells1Red,1:cells(2),1:cells3) = sum(vectorField_fourier(1:3,1:cells1Red,1:cells(2),1:cells3) &
*conjg(-xi1st),1)
end subroutine utilities_fourierVectorDivergence
@ -764,11 +801,12 @@ subroutine utilities_fourierVectorGradient()
integer :: i, j, k, m, n
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells1Red
do m = 1, 3; do n = 1, 3
tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k)
enddo; enddo
enddo; enddo; enddo
end do; end do
end do; end do; end do
end subroutine utilities_fourierVectorGradient
@ -780,9 +818,10 @@ subroutine utilities_fourierTensorDivergence()
integer :: i, j, k
do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells1Red
vectorField_fourier(:,i,j,k) = matmul(tensorField_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k)))
enddo; enddo; enddo
end do; end do; end do
end subroutine utilities_fourierTensorDivergence
@ -812,9 +851,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field
call homogenization_mechanical_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) ! calculate P field
call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3) ! calculate P field
if (.not. terminallyIll) &
call homogenization_thermal_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3)
if (.not. terminallyIll) &
call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
@ -884,11 +923,10 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_calculateRate
if (heterogeneous) then
utilities_calculateRate = (field-field0) / dt
else
utilities_calculateRate = spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3)
endif
utilities_calculateRate = merge((field-field0) / dt, &
spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3), &
heterogeneous)
end function utilities_calculateRate
@ -980,6 +1018,7 @@ end function utilities_getFreqDerivative
subroutine utilities_updateCoords(F)
real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: F
real(pReal), dimension(3, cells(1),cells(2),cells3) :: IPcoords
real(pReal), dimension(3, cells(1),cells(2),cells3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI)
real(pReal), dimension(3, cells(1)+1,cells(2)+1,cells3+1) :: nodeCoords
@ -1010,20 +1049,23 @@ subroutine utilities_updateCoords(F)
1, 1, 1, &
0, 1, 1 ], [3,8])
step = geomSize/real(cells, pReal)
!--------------------------------------------------------------------------------------------------
! integration in Fourier space to get fluctuations of cell center discplacements
tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F
call utilities_FFTtensorForward()
do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red
!$OMP PARALLEL DO
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells1Red
if (any([i,j,k+cells3Offset] /= 1)) then
vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) &
/ sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal)
else
vectorField_fourier(1:3,i,j,k) = cmplx(0.0,0.0,pReal)
endif
enddo; enddo; enddo
end if
end do; end do; end do
!$OMP END PARALLEL DO
call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real)
@ -1041,7 +1083,7 @@ subroutine utilities_updateCoords(F)
rank_b = modulo(worldrank-1_MPI_INTEGER_KIND,worldsize)
! send bottom layer to process below
call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(1),err_MPI)
call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(1),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Irecv(IPfluct_padded(:,:,:,cells3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@ -1049,7 +1091,7 @@ subroutine utilities_updateCoords(F)
! send top layer to process above
call MPI_Isend(IPfluct_padded(:,:,:,cells3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(4),err_MPI)
call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(4),err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Waitall(4,request,status,err_MPI)

106
src/homogenization.f90
View File

@ -222,14 +222,13 @@ end subroutine homogenization_init
!--------------------------------------------------------------------------------------------------
!> @brief
!--------------------------------------------------------------------------------------------------
subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
real(pReal), intent(in) :: Delta_t !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer, intent(in) :: &
cell_start, cell_end
integer :: &
NiterationMPstate, &
ip, & !< integration point number
el, & !< element number
co, ce, ho, en
logical :: &
converged
@ -237,45 +236,42 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving
doneAndHappy
!$OMP PARALLEL DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
!$OMP PARALLEL DO PRIVATE(en,ho,co,NiterationMPstate,converged,doneAndHappy)
do ce = cell_start, cell_end
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip
en = material_homogenizationEntry(ce)
ho = material_homogenizationID(ce)
en = material_homogenizationEntry(ce)
ho = material_homogenizationID(ce)
call phase_restore(ce,.false.) ! wrong name (is more a forward function)
call phase_restore(ce,.false.) ! wrong name (is more a forward function)
if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en)
call damage_partition(ce)
if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en)
call damage_partition(ce)
doneAndHappy = [.false.,.true.]
doneAndHappy = [.false.,.true.]
NiterationMPstate = 0
convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) &
.and. NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1
NiterationMPstate = 0
convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) &
.and. NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1
call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
if (converged) then
doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
converged = all(doneAndHappy)
else
doneAndHappy = [.true.,.false.]
endif
enddo convergenceLooping
call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))])
if (converged) then
doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
converged = all(doneAndHappy)
else
doneAndHappy = [.true.,.false.]
end if
end do convergenceLooping
converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))])
if (.not. converged) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
endif
enddo
enddo
if (.not. converged) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
end if
end do
!$OMP END PARALLEL DO
end subroutine homogenization_mechanical_response
@ -284,31 +280,27 @@ end subroutine homogenization_mechanical_response
!--------------------------------------------------------------------------------------------------
!> @brief
!--------------------------------------------------------------------------------------------------
subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end)
real(pReal), intent(in) :: Delta_t !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer, intent(in) :: &
cell_start, cell_end
integer :: &
ip, & !< integration point number
el, & !< element number
co, ce, ho
!$OMP PARALLEL DO PRIVATE(ho,ce)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
!$OMP PARALLEL DO PRIVATE(ho)
do ce = cell_start, cell_end
if (terminallyIll) continue
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip
ho = material_homogenizationID(ce)
call thermal_partition(ce)
do co = 1, homogenization_Nconstituents(ho)
if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
endif
enddo
enddo
enddo
ho = material_homogenizationID(ce)
call thermal_partition(ce)
do co = 1, homogenization_Nconstituents(ho)
if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
end if
end do
end do
!$OMP END PARALLEL DO
end subroutine homogenization_thermal_response
@ -331,11 +323,11 @@ subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolvin
elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip
ho = material_homogenizationID(ce)
do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el)
ho = material_homogenizationID(ce)
do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el)
enddo
call mechanical_homogenize(Delta_t,ce)
call mechanical_homogenize(Delta_t,ce)
enddo IpLooping3
enddo elementLooping3
!$OMP END PARALLEL DO

13
src/homogenization_mechanical.f90
View File

@ -131,20 +131,15 @@ module subroutine mechanical_homogenize(Delta_t,ce)
integer :: co
homogenization_P(1:3,1:3,ce) = phase_P(1,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce)
homogenization_P(1:3,1:3,ce) = phase_P(1,ce)*material_v(1,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce)*material_v(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) &
+ phase_P(co,ce)
+ phase_P(co,ce)*material_v(co,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
+ phase_mechanical_dPdF(Delta_t,co,ce)
+ phase_mechanical_dPdF(Delta_t,co,ce)*material_v(co,ce)
end do
homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) &
/ real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
/ real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end subroutine mechanical_homogenize

18
src/homogenization_thermal.f90
View File

@ -105,13 +105,11 @@ module function homogenization_mu_T(ce) result(mu)
integer :: co
mu = phase_mu_T(1,ce)
mu = phase_mu_T(1,ce)*material_v(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
mu = mu + phase_mu_T(co,ce)
mu = mu + phase_mu_T(co,ce)*material_v(co,ce)
end do
mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function homogenization_mu_T
@ -126,13 +124,11 @@ module function homogenization_K_T(ce) result(K)
integer :: co
K = phase_K_T(1,ce)
K = phase_K_T(1,ce)*material_v(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
K = K + phase_K_T(co,ce)
K = K + phase_K_T(co,ce)*material_v(co,ce)
end do
K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function homogenization_K_T
@ -147,13 +143,11 @@ module function homogenization_f_T(ce) result(f)
integer :: co
f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))*material_v(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))*material_v(co,ce)
end do
f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function homogenization_f_T

92
src/material.f90
View File

@ -17,16 +17,17 @@ module material
implicit none
private
type :: tRotationContainer
type(tRotation), dimension(:), allocatable :: data
end type
type :: tTensorContainer
type, public :: tRotationContainer
type(tRotation), dimension(:), allocatable :: data
end type tRotationContainer
type, public :: tTensorContainer
real(pReal), dimension(:,:,:), allocatable :: data
end type
end type tTensorContainer
type(tRotationContainer), dimension(:), allocatable :: material_O_0
type(tTensorContainer), dimension(:), allocatable :: material_F_i_0
type(tRotationContainer), dimension(:), allocatable, public, protected :: material_O_0
type(tTensorContainer), dimension(:), allocatable, public, protected :: material_F_i_0
integer, dimension(:), allocatable, public, protected :: &
homogenization_Nconstituents !< number of grains in each homogenization
@ -37,20 +38,17 @@ module material
material_name_phase, & !< name of each phase
material_name_homogenization !< name of each homogenization
integer, dimension(:), allocatable, public, protected :: & ! (elem)
material_homogenizationID, & !< per cell TODO: material_ID_homogenization
material_homogenizationEntry !< per cell TODO: material_entry_homogenization
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt, & !< phase ID of each element TODO: remove
material_phaseID, & !< per (constituent,cell) TODO: material_ID_phase
material_phaseEntry !< per (constituent,cell) TODO: material_entry_phase
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem)
material_phaseMemberAt !TODO: remove
integer, dimension(:), allocatable, public, protected :: & ! (cell)
material_homogenizationID, & ! TODO: rename to material_ID_homogenization
material_homogenizationEntry ! TODO: rename to material_entry_homogenization
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,cell)
material_phaseID, & ! TODO: rename to material_ID_phase
material_phaseEntry ! TODO: rename to material_entry_phase
real(pReal), dimension(:,:), allocatable, public, protected :: &
material_v ! fraction
public :: &
tTensorContainer, &
tRotationContainer, &
material_F_i_0, &
material_O_0, &
material_init
contains
@ -97,11 +95,12 @@ subroutine parse()
counterPhase, &
counterHomogenization
real(pReal) :: &
frac
real(pReal) :: v
integer :: &
el, ip, co, ma, &
h, ce
el, ip, &
ho, ph, &
co, ce, &
ma
materials => config_material%get('material')
phases => config_material%get('phase')
@ -118,51 +117,52 @@ subroutine parse()
#endif
allocate(homogenization_Nconstituents(homogenizations%length))
do h=1, homogenizations%length
homogenization => homogenizations%get(h)
homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents')
do ho=1, homogenizations%length
homogenization => homogenizations%get(ho)
homogenization_Nconstituents(ho) = homogenization%get_asInt('N_constituents')
end do
homogenization_maxNconstituents = maxval(homogenization_Nconstituents)
allocate(counterPhase(phases%length),source=0)
allocate(counterHomogenization(homogenizations%length),source=0)
allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0)
allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0)
allocate(material_homogenizationID(discretization_Ncells),source=0)
allocate(material_homogenizationEntry(discretization_Ncells),source=0)
allocate(material_phaseID(homogenization_maxNconstituents,discretization_Ncells),source=0)
allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_Ncells),source=0)
allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
allocate(material_v(homogenization_maxNconstituents,discretization_Ncells),source=0.0_pReal)
do el = 1, discretization_Nelems
material => materials%get(discretization_materialAt(el))
constituents => material%get('constituents')
material => materials%get(discretization_materialAt(el))
ho = homogenizations%getIndex(material%get_asString('homogenization'))
do ip = 1, discretization_nIPs
ce = (el-1)*discretization_nIPs + ip
material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization'))
counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1
material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce))
material_homogenizationID(ce) = ho
counterHomogenization(ho) = counterHomogenization(ho) + 1
material_homogenizationEntry(ce) = counterHomogenization(ho)
end do
frac = 0.0_pReal
constituents => material%get('constituents')
do co = 1, constituents%length
constituent => constituents%get(co)
frac = frac + constituent%get_asFloat('v')
material_phaseAt(co,el) = phases%getIndex(constituent%get_asString('phase'))
v = constituent%get_asFloat('v')
ph = phases%getIndex(constituent%get_asString('phase'))
do ip = 1, discretization_nIPs
ce = (el-1)*discretization_nIPs + ip
counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1
material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el))
material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el))
material_phaseID(co,ce) = material_phaseAt(co,el)
material_phaseID(co,ce) = ph
counterPhase(ph) = counterPhase(ph) + 1
material_phaseEntry(co,ce) = counterPhase(ph)
material_v(co,ce) = v
end do
end do
if (dNeq(frac,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent')
if (dNeq(sum(material_v(1:constituents%length,ce)),1.0_pReal,1.e-9_pReal)) &
call IO_error(153,ext_msg='constituent')
end do

31
src/math.f90
View File

@ -262,9 +262,8 @@ pure function math_identity4th()
math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
enddo
#else
do i=1,3; do j=1,3; do k=1,3; do l=1,3
forall(i=1:3, j=1:3, k=1:3, l=1:3) &
math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
enddo; enddo; enddo; enddo
#endif
end function math_identity4th
@ -338,9 +337,7 @@ pure function math_outer(A,B)
math_outer(i,j) = A(i)*B(j)
enddo
#else
do i=1,size(A,1); do j=1,size(B,1)
math_outer(i,j) = A(i)*B(j)
enddo; enddo
forall(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j)
#endif
end function math_outer
@ -387,9 +384,7 @@ pure function math_mul3333xx33(A,B)
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
enddo
#else
do i=1,3; do j=1,3
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
enddo; enddo
forall (i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
#endif
end function math_mul3333xx33
@ -411,9 +406,7 @@ pure function math_mul3333xx3333(A,B)
math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
enddo
#else
do i=1,3; do j=1,3; do k=1,3; do l=1,3
math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
enddo; enddo; enddo; enddo
forall(i=1:3, j=1:3, k=1:3, l=1:3) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
#endif
end function math_mul3333xx3333
@ -752,9 +745,7 @@ pure function math_3333to99(m3333)
math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
enddo
#else
do i=1,9; do j=1,9
math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
enddo; enddo
forall(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
#endif
end function math_3333to99
@ -775,9 +766,7 @@ pure function math_99to3333(m99)
math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
enddo
#else
do i=1,9; do j=1,9
math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
enddo; enddo
forall(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
#endif
end function math_99to3333
@ -810,9 +799,7 @@ pure function math_sym3333to66(m3333,weighted)
math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
enddo
#else
do i=1,6; do j=1,6
math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
enddo; enddo
forall(i=1:6, j=1:6) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
#endif
end function math_sym3333to66
@ -950,9 +937,7 @@ pure function math_3333toVoigt66_stiffness(C) result(C_tilde)
C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
end do
#else
do i=1,6; do j=1,6
C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
end do; end do
forall(i=1:6, j=1:6) C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
#endif
end function math_3333toVoigt66_stiffness

2
src/mesh/FEM_utilities.f90
View File

@ -150,7 +150,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
print'(/,1x,a)', '... evaluating constitutive response ......................................'
call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field
call homogenization_mechanical_response(timeinc,1,mesh_maxNips*mesh_NcpElems) ! calculate P field
if (.not. terminallyIll) &
call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])
cutBack = .false.

37
src/phase.f90
View File

@ -228,15 +228,15 @@ module phase
end function phase_thermal_constitutive
module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
module function phase_damage_constitutive(Delta_t,co,ce) result(converged_)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co, ip, el
integer, intent(in) :: co, ce
logical :: converged_
end function phase_damage_constitutive
module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co, ip, el
integer, intent(in) :: co, ce
logical :: converged_
end function phase_mechanical_constitutive
@ -264,19 +264,18 @@ module phase
real(pReal) :: f
end function phase_f_T
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el)
integer, intent(in) :: &
ph, &
i, &
e
ip, &
el
type(tRotationContainer), dimension(:), intent(in) :: orientation
end subroutine plastic_nonlocal_updateCompatibility
module subroutine plastic_dependentState(co,ip,el)
module subroutine plastic_dependentState(ph,en)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
ph, &
en
end subroutine plastic_dependentState
@ -437,6 +436,8 @@ subroutine phase_allocateState(state, &
allocate(state%dotState (sizeDotState,NEntries), source=0.0_pReal)
allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal)
state%deltaState2 => state%state(state%offsetDeltaState+1: &
state%offsetDeltaState+state%sizeDeltaState,:)
end subroutine phase_allocateState
@ -501,8 +502,8 @@ subroutine crystallite_init()
ce, &
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el !< counter in element loop
el, & !< counter in element loop
en, ph
class(tNode), pointer :: &
num_crystallite, &
phases
@ -541,13 +542,15 @@ subroutine crystallite_init()
phases => config_material%get('phase')
!$OMP PARALLEL DO PRIVATE(ce)
!$OMP PARALLEL DO PRIVATE(ce,ph,en)
do el = 1, discretization_Nelems
do ip = 1, discretization_nIPs
ce = (el-1)*discretization_nIPs + ip
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
en = material_phaseEntry(co,ce)
ph = material_phaseID(co,ce)
call crystallite_orientations(co,ip,el)
call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
call plastic_dependentState(ph,en) ! update dependent state variables to be consistent with basic states
end do
end do
end do
@ -575,8 +578,8 @@ subroutine crystallite_orientations(co,ip,el)
call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en))))
if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(phase_O,material_phaseAt(1,el),ip,el)
if (plasticState(material_phaseID(1,(el-1)*discretization_nIPs + ip))%nonlocal) &
call plastic_nonlocal_updateCompatibility(phase_O,material_phaseID(1,(el-1)*discretization_nIPs + ip),ip,el)
end subroutine crystallite_orientations

9
src/phase_damage.f90
View File

@ -127,21 +127,20 @@ end subroutine damage_init
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
module function phase_damage_constitutive(Delta_t,co,ce) result(converged_)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
co, &
ip, &
el
ce
logical :: converged_
integer :: &
ph, en
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
ph = material_phaseID(co,ce)
en = material_phaseEntry(co,ce)
converged_ = .not. integrateDamageState(Delta_t,ph,en)

130
src/phase_mechanical.f90
View File

@ -79,11 +79,8 @@ submodule(phase) mechanical
en
end subroutine plastic_isotropic_LiAndItsTangent
module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
module function plastic_dotState(subdt,ph,en) result(dotState)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
en
real(pReal), intent(in) :: &
@ -365,13 +362,11 @@ end subroutine mechanical_results
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction
!--------------------------------------------------------------------------------------------------
function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
real(pReal), dimension(3,3), intent(in) :: F,subFp0,subFi0
real(pReal), intent(in) :: Delta_t
integer, intent(in):: el, & ! element index
ip, & ! integration point index
co ! grain index
integer, intent(in) :: ph, en
real(pReal), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep
invFp_new, & ! inverse of Fp_new
@ -419,19 +414,13 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
ierr, & ! error indicator for LAPACK
o, &
p, &
ph, &
en, &
jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update
logical :: error,broken
broken = .true.
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
call plastic_dependentState(co,ip,el)
call plastic_dependentState(ph,en)
Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess
Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess
@ -579,37 +568,32 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
ph, &
en
logical :: &
broken
integer :: &
NiterationState, & !< number of iterations in state loop
ph, &
en, &
sizeDotState
real(pReal) :: &
zeta
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
r, & ! state residuum
dotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState,2) :: &
dotState_last
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
@ -620,10 +604,10 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1)
dotState_last(1:sizeDotState,1) = dotState
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
if(broken) exit iteration
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) exit iteration
zeta = damper(dotState,dotState_last(1:sizeDotState,1),dotState_last(1:sizeDotState,2))
@ -672,31 +656,26 @@ end function integrateStateFPI
!--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
ph, &
en !< grain index in grain loop
logical :: &
broken
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
integer :: &
ph, &
en, &
sizeDotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
@ -706,7 +685,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
broken = plastic_deltaState(ph,en)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
end function integrateStateEuler
@ -714,32 +693,27 @@ end function integrateStateEuler
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
ph, &
en
logical :: &
broken
integer :: &
ph, &
en, &
sizeDotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
r, &
dotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
@ -751,10 +725,10 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
broken = plastic_deltaState(ph,en)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
broken = .not. converged(r + 0.5_pReal * dotState * Delta_t, &
@ -767,12 +741,12 @@ end function integrateStateAdaptiveEuler
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the classic Runge Kutta method
!---------------------------------------------------------------------------------------------------
function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co,ip,el
integer, intent(in) :: ph, en
logical :: broken
real(pReal), dimension(3,3), parameter :: &
@ -787,7 +761,7 @@ function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal]
broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C)
broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C)
end function integrateStateRK4
@ -795,12 +769,12 @@ end function integrateStateRK4
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method
!---------------------------------------------------------------------------------------------------
function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co,ip,el
integer, intent(in) :: ph, en
logical :: broken
real(pReal), dimension(5,5), parameter :: &
@ -822,7 +796,7 @@ function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) re
13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal]
broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB)
broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
end function integrateStateRKCK45
@ -831,7 +805,7 @@ end function integrateStateRKCK45
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method
!--------------------------------------------------------------------------------------------------
function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB) result(broken)
function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
@ -840,28 +814,23 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
real(pReal), dimension(:), intent(in) :: B, C
real(pReal), dimension(:), intent(in), optional :: DB
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
ph, &
en
logical :: broken
integer :: &
stage, & ! stage index in integration stage loop
n, &
ph, &
en, &
sizeDotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState,size(B)) :: &
plastic_RKdotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
@ -879,10 +848,10 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ dotState * Delta_t
broken = integrateStress(F_0 + (F-F_0) * Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage),co,ip,el)
broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage), ph,en)
if(broken) exit
dotState = plastic_dotState(Delta_t*C(stage), co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t*C(stage), ph,en)
if (any(IEEE_is_NaN(dotState))) exit
enddo
@ -904,7 +873,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
broken = plastic_deltaState(ph,en)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
end function integrateStateRK
@ -1006,13 +975,12 @@ end subroutine mechanical_forward
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
co, &
ip, &
el
ce
logical :: converged_
real(pReal) :: &
@ -1028,16 +996,15 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
subLi0, &
subF0, &
subF
real(pReal), dimension(:), allocatable :: subState0
real(pReal), dimension(plasticState(material_phaseID(co,ce))%sizeState) :: subState0
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
sizeDotState = plasticState(ph)%sizeDotState
ph = material_phaseID(co,ce)
en = material_phaseEntry(co,ce)
subState0 = plasticState(ph)%state0(:,en)
subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
allocate(subState0,source=plasticState(ph)%State0(:,en))
subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en)
@ -1070,7 +1037,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
subStep = num%subStepSizeCryst * subStep
phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0
phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0
phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use?
phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)
if (subStep < 1.0_pReal) then ! actual (not initial) cutback
phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0
phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0
@ -1082,9 +1049,10 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
!--------------------------------------------------------------------------------------------------
! prepare for integration
if (todo) then
sizeDotState = plasticState(ph)%sizeDotState
subF = subF0 &
+ subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el)
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en)
endif
enddo cutbackLooping

2
src/phase_mechanical_elastic.f90
View File

@ -48,7 +48,7 @@ module subroutine elastic_init(phases)
prm%C_11 = polynomial(elastic%asDict(),'C_11','T')
prm%C_12 = polynomial(elastic%asDict(),'C_12','T')
prm%C_44 = polynomial(elastic%asDict(),'C_44','T')
if (any(phase_lattice(ph) == ['hP','tI'])) then
prm%C_13 = polynomial(elastic%asDict(),'C_13','T')
prm%C_33 = polynomial(elastic%asDict(),'C_33','T')

78
src/phase_mechanical_plastic.f90
View File

@ -110,29 +110,35 @@ submodule(phase:mechanical) plastic
end subroutine nonlocal_LpAndItsTangent
module subroutine isotropic_dotState(Mp,ph,en)
module function isotropic_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine isotropic_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function isotropic_dotState
module subroutine phenopowerlaw_dotState(Mp,ph,en)
module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine phenopowerlaw_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function phenopowerlaw_dotState
module subroutine plastic_kinehardening_dotState(Mp,ph,en)
module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine plastic_kinehardening_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function plastic_kinehardening_dotState
module subroutine dislotwin_dotState(Mp,T,ph,en)
real(pReal), dimension(3,3), intent(in) :: &
@ -144,26 +150,24 @@ submodule(phase:mechanical) plastic
en
end subroutine dislotwin_dotState
module subroutine dislotungsten_dotState(Mp,T,ph,en)
module function dislotungsten_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
ph, &
en
end subroutine dislotungsten_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function dislotungsten_dotState
module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el)
module subroutine nonlocal_dotState(Mp,timestep,ph,en)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
en, &
ip, & !< current integration point
el !< current element number
en
end subroutine nonlocal_dotState
module subroutine dislotwin_dependentState(T,ph,en)
@ -180,12 +184,10 @@ submodule(phase:mechanical) plastic
en
end subroutine dislotungsten_dependentState
module subroutine nonlocal_dependentState(ph, en, ip, el)
module subroutine nonlocal_dependentState(ph,en)
integer, intent(in) :: &
ph, &
en, &
ip, & !< current integration point
el !< current element number
en
end subroutine nonlocal_dependentState
module subroutine plastic_kinehardening_deltaState(Mp,ph,en)
@ -295,12 +297,9 @@ end subroutine plastic_LpAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
module function plastic_dotState(subdt,ph,en) result(dotState)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
en
real(pReal), intent(in) :: &
@ -318,48 +317,41 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
plasticType: select case (phase_plasticity(ph))
case (PLASTIC_ISOTROPIC_ID) plasticType
call isotropic_dotState(Mp,ph,en)
dotState = isotropic_dotState(Mp,ph,en)
case (PLASTIC_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_dotState(Mp,ph,en)
dotState = phenopowerlaw_dotState(Mp,ph,en)
case (PLASTIC_KINEHARDENING_ID) plasticType
call plastic_kinehardening_dotState(Mp,ph,en)
dotState = plastic_kinehardening_dotState(Mp,ph,en)
case (PLASTIC_DISLOTWIN_ID) plasticType
call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en)
dotState = plasticState(ph)%dotState(:,en)
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en)
dotState = dislotungsten_dotState(Mp,ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType
call nonlocal_dotState(Mp,subdt,ph,en,ip,el)
call nonlocal_dotState(Mp,subdt,ph,en)
dotState = plasticState(ph)%dotState(:,en)
end select plasticType
end if
dotState = plasticState(ph)%dotState(:,en)
end function plastic_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different plasticity constitutive models
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dependentState(co, ip, el)
module subroutine plastic_dependentState(ph,en)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: &
ph, &
en
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
plasticType: select case (phase_plasticity(ph))
case (PLASTIC_DISLOTWIN_ID) plasticType
@ -369,7 +361,7 @@ module subroutine plastic_dependentState(co, ip, el)
call dislotungsten_dependentState(ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType
call nonlocal_dependentState(ph,en,ip,el)
call nonlocal_dependentState(ph,en)
end select plasticType
@ -387,10 +379,9 @@ module function plastic_deltaState(ph, en) result(broken)
en
logical :: broken
real(pReal), dimension(3,3) :: &
real(pReal), dimension(3,3) :: &
Mp
integer :: &
myOffset, &
mySize
@ -415,10 +406,9 @@ module function plastic_deltaState(ph, en) result(broken)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
if (.not. broken) then
myOffset = plasticState(ph)%offsetDeltaState
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = &
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en)
plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) &
+ plasticState(ph)%deltaState(1:mySize,en)
end if
end select

57
src/phase_mechanical_plastic_dislotungsten.f90
View File

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

44
src/phase_mechanical_plastic_isotropic.f90
View File

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

56
src/phase_mechanical_plastic_kinehardening.f90
View File

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

152
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -15,6 +15,9 @@ submodule(phase:plastic) nonlocal
type :: tGeometry
real(pReal), dimension(:), allocatable :: V_0
integer, dimension(:,:,:), allocatable :: IPneighborhood
real(pReal), dimension(:,:), allocatable :: IParea, IPcoordinates
real(pReal), dimension(:,:,:), allocatable :: IPareaNormal
end type tGeometry
type(tGeometry), dimension(:), allocatable :: geom
@ -123,8 +126,10 @@ submodule(phase:plastic) nonlocal
type :: tNonlocalDependentState
real(pReal), allocatable, dimension(:,:) :: &
tau_pass, &
tau_Back
tau_pass, &
tau_Back
real(pReal), allocatable, dimension(:,:,:,:,:) :: &
compatibility
end type tNonlocalDependentState
type :: tNonlocalState
@ -160,7 +165,6 @@ submodule(phase:plastic) nonlocal
state0
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters
type(tNonlocalDependentState), dimension(:), allocatable :: dependentState
contains
@ -227,8 +231,8 @@ module function plastic_nonlocal_init() result(myPlasticity)
st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
plasticState(ph)%nonlocal = pl%get_asBool('flux',defaultVal=.True.)
#if defined (__GFORTRAN__)
@ -406,6 +410,10 @@ module function plastic_nonlocal_init() result(myPlasticity)
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention
allocate(geom(ph)%V_0(Nmembers))
allocate(geom(ph)%IPneighborhood(3,nIPneighbors,Nmembers))
allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers))
allocate(geom(ph)%IParea(nIPneighbors,Nmembers))
allocate(geom(ph)%IPcoordinates(3,Nmembers))
call storeGeometry(ph)
if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
@ -489,6 +497,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,Nmembers),source=0.0_pReal)
end associate
if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
@ -543,13 +552,11 @@ end function plastic_nonlocal_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure
!--------------------------------------------------------------------------------------------------
module subroutine nonlocal_dependentState(ph, en, ip, el)
module subroutine nonlocal_dependentState(ph, en)
integer, intent(in) :: &
ph, &
en, &
ip, &
el
en
integer :: &
no, & !< neighbor offset
@ -650,12 +657,12 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
nRealNeighbors = 0.0_pReal
neighbor_rhoTotal = 0.0_pReal
do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
if (neighbor_el > 0 .and. neighbor_ip > 0) then
if (material_phaseAt(1,neighbor_el) == ph) then
neighbor_el = geom(ph)%IPneighborhood(1,n,en)
neighbor_ip = geom(ph)%IPneighborhood(2,n,en)
if (neighbor_el > 0 .and. neighbor_ip > 0) then
if (material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) == ph) then
no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip)
nRealNeighbors = nRealNeighbors + 1.0_pReal
rho_neighbor0 = getRho0(ph,no)
@ -665,11 +672,11 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2)
neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2)
connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) &
- discretization_IPcoords(1:3,el+neighbor_ip-1))
normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el))
connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%IPcoordinates(1:3,no) &
- geom(ph)%IPcoordinates(1:3,en))
normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(1:3,n,en))
if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image
connection_latticeConf(1:3,n) = normal_latticeConf * IPvolume(ip,el)/IParea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/geom(ph)%IParea(n,en) ! instead take the surface normal scaled with the diameter of the cell
else
! local neighbor or different lattice structure or different constitution instance -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pReal
@ -939,7 +946,7 @@ end subroutine plastic_nonlocal_deltaState
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
module subroutine nonlocal_dotState(Mp,timestep, &
ph,en,ip,el)
ph,en)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
@ -947,9 +954,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
en, &
ip, & !< current integration point
el !< current element number
en
integer :: &
c, & !< character of dislocation
@ -1099,7 +1104,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(timestep, ph,en,ip,el) &
rhoDot = rhoDotFlux(timestep, ph,en) &
+ rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation &
@ -1110,7 +1115,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
.or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at ph ',ph,' en ',en
print'(a)', '<< CONST >> enforcing cutback !!!'
end if
#endif
@ -1129,20 +1134,17 @@ end subroutine nonlocal_dotState
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
#if __INTEL_COMPILER >= 2020
non_recursive function rhoDotFlux(timestep,ph,en,ip,el)
non_recursive function rhoDotFlux(timestep,ph,en)
#else
function rhoDotFlux(timestep,ph,en,ip,el)
function rhoDotFlux(timestep,ph,en)
#endif
real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
en, &
ip, & !< current integration point
el !< current element number
en
integer :: &
neighbor_ph, & !< phase of my neighbor's plasticity
ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation
n, & !< index of my current neighbor
@ -1215,14 +1217,14 @@ function rhoDotFlux(timestep,ph,en,ip,el)
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
.and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
> geom(ph)%V_0(en)/ maxval(geom(ph)%IParea(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at ph ',ph,' en ',en
print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
.and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
> geom(ph)%V_0(en) / maxval(geom(ph)%IParea(:,en))), &
' at a timestep of ',timestep
print*, '<< CONST >> enforcing cutback !!!'
end if
@ -1245,19 +1247,18 @@ function rhoDotFlux(timestep,ph,en,ip,el)
neighbors: do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
neighbor_n = IPneighborhood(3,n,ip,el)
np = material_phaseAt(1,neighbor_el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
neighbor_el = geom(ph)%IPneighborhood(1,n,en)
neighbor_ip = geom(ph)%IPneighborhood(2,n,en)
neighbor_n = geom(ph)%IPneighborhood(3,n,en)
np = material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip)
no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip)
opposite_neighbor = n + mod(n,2) - mod(n+1,2)
opposite_el = IPneighborhood(1,opposite_neighbor,ip,el)
opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el)
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el)
opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en)
opposite_ip = geom(ph)%IPneighborhood(2,opposite_neighbor,en)
opposite_n = geom(ph)%IPneighborhood(3,opposite_neighbor,en)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_ph = material_phaseAt(1,neighbor_el)
neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no)
neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F)
@ -1276,12 +1277,12 @@ function rhoDotFlux(timestep,ph,en,ip,el)
!* The entering flux from my neighbor will be distributed on my slip systems according to the
!* compatibility
if (neighbor_n > 0) then
if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTIC_NONLOCAL_ID .and. &
any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then
if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID .and. &
any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pReal)) then
forall (s = 1:ns, t = 1:4)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_ph),no)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_ph),no),0.0_pReal)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,np),no)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,np),no),0.0_pReal)
endforall
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min &
@ -1301,12 +1302,12 @@ function rhoDotFlux(timestep,ph,en,ip,el)
.and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density
lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) &
* math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (compatibility(c,:,s,n,ip,el) > 0.0_pReal) &
where (dependentState(ph)%compatibility(c,:,s,n,en) > 0.0_pReal) &
rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) &
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to equally signed mobile dislocation type
where (compatibility(c,:,s,n,ip,el) < 0.0_pReal) &
+ lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to equally signed mobile dislocation type
where (dependentState(ph)%compatibility(c,:,s,n,en) < 0.0_pReal) &
rhoDotFlux(:,topp) = rhoDotFlux(:,topp) &
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to opposite signed mobile dislocation type
+ lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to opposite signed mobile dislocation type
end if
end do
@ -1322,28 +1323,28 @@ function rhoDotFlux(timestep,ph,en,ip,el)
!* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
!* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
if (opposite_n > 0) then
if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then
if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID) then
normal_me2neighbor_defConf = math_det33(Favg) &
* matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
* matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) &
/ math_det33(my_Fe) ! interface normal in my lattice configuration
area = IParea(n,ip,el) * norm2(normal_me2neighbor)
area = geom(ph)%IParea(n,en) * norm2(normal_me2neighbor)
normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length
do s = 1,ns
do t = 1,4
c = (t + 1) / 2
if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density
transmissivity = sum(compatibility(c,:,s,n,ip,el)**2) ! overall transmissivity from this slip system to my neighbor
transmissivity = sum(dependentState(ph)%compatibility(c,:,s,n,en)**2) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal
end if
lineLength = my_rhoSgl0(s,t) * v0(s,t) &
* math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / geom(ph)%V_0(en) ! subtract dislocation flux from current type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
+ lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) &
+ lineLength / geom(ph)%V_0(en) * (1.0_pReal - transmissivity) &
* sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
end if
end do
@ -1364,14 +1365,14 @@ end function rhoDotFlux
! plane normals and signed cosine of the angle between the slip directions. Only the largest values
! that sum up to a total of 1 are considered, all others are set to zero.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el)
type(tRotationContainer), dimension(:), intent(in) :: &
orientation ! crystal orientation
integer, intent(in) :: &
ph, &
i, &
e
ip, &
el
integer :: &
n, & ! neighbor index
@ -1397,16 +1398,16 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
associate(prm => param(ph))
ns = prm%sum_N_sl
en = material_phaseMemberAt(1,i,e)
en = material_phaseEntry(1,(el-1)*discretization_nIPs + ip)
!*** start out fully compatible
my_compatibility = 0.0_pReal
forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
neighbors: do n = 1,nIPneighbors
neighbor_e = IPneighborhood(1,n,i,e)
neighbor_i = IPneighborhood(2,n,i,e)
neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
neighbor_phase = material_phaseAt(1,neighbor_e)
neighbor_e = IPneighborhood(1,n,ip,el)
neighbor_i = IPneighborhood(2,n,ip,el)
neighbor_me = material_phaseEntry(1,(neighbor_e-1)*discretization_nIPs + neighbor_i)
neighbor_phase = material_phaseID(1,(neighbor_e-1)*discretization_nIPs + neighbor_i)
if (neighbor_e <= 0 .or. neighbor_i <= 0) then
!* FREE SURFACE
@ -1465,7 +1466,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
end do neighbors
compatibility(:,:,:,:,i,e) = my_compatibility
dependentState(ph)%compatibility(:,:,:,:,material_phaseEntry(1,(el-1)*discretization_nIPs + ip)) = my_compatibility
end associate
@ -1756,14 +1757,29 @@ subroutine storeGeometry(ph)
integer, intent(in) :: ph
integer :: ce, co
integer :: ce, co, nCell
real(pReal), dimension(:), allocatable :: V
integer, dimension(:,:,:), allocatable :: neighborhood
real(pReal), dimension(:,:), allocatable :: area, coords
real(pReal), dimension(:,:,:), allocatable :: areaNormal
nCell = product(shape(IPvolume))
V = reshape(IPvolume,[nCell])
neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell])
area = reshape(IParea,[nIPneighbors,nCell])
areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell])
coords = reshape(discretization_IPcoords,[3,nCell])
V = reshape(IPvolume,[product(shape(IPvolume))])
do ce = 1, size(material_homogenizationEntry,1)
do co = 1, homogenization_maxNconstituents
if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
if (material_phaseID(co,ce) == ph) then
geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
geom(ph)%IPneighborhood(:,:,material_phaseEntry(co,ce)) = neighborhood(:,:,ce)
geom(ph)%IParea(:,material_phaseEntry(co,ce)) = area(:,ce)
geom(ph)%IPareaNormal(:,:,material_phaseEntry(co,ce)) = areaNormal(:,:,ce)
geom(ph)%IPcoordinates(:,material_phaseEntry(co,ce)) = coords(:,ce)
end if
end do
end do

70
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

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

21
src/phase_thermal.f90
View File

@ -96,17 +96,24 @@ module subroutine thermal_init(phases)
do ph = 1, phases%length
Nmembers = count(material_phaseID == ph)
allocate(current(ph)%T(Nmembers),source=300.0_pReal)
allocate(current(ph)%T(Nmembers),source=T_ROOM)
allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
phase => phases%get(ph)
thermal => phase%get('thermal',defaultVal=emptyDict)
param(ph)%C_p = thermal%get_asFloat('C_p',defaultVal=0.0_pReal) ! ToDo: make mandatory?
param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal) ! ToDo: make mandatory?
param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal) ! ToDo: depends on symmtery
param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length
! ToDo: temperature dependency of K and C_p
if (thermal%length > 0) then
param(ph)%C_p = thermal%get_asFloat('C_p')
param(ph)%K(1,1) = thermal%get_asFloat('K_11')
if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asFloat('K_33')
param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length
else
thermal_Nsources(ph) = 0
end if
allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
enddo

2
src/prec.f90
View File

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