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

Merge branch 'Fortran-cleaning' into 'development' 

Clarifying structure

See merge request damask/DAMASK!368
This commit is contained in:
Sharan Roongta 2021-04-13 11:08:18 +00:00
parent 90e33b14bf 887524bcc1
commit 1be1a72a09
28 changed files with 905 additions and 1054 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 1298124143e7e2901d0b9c2e79ab6388cb78a1e3 Subproject commit 1490f97417664d6ae11d7ceafb6b98c9cc2dade1

3
examples/.gitignore vendored
View File

@ -3,3 +3,6 @@
*.xdmf *.xdmf
*.sta *.sta
*.vt* *.vt*
*.out
*.sts
*.t16

429
examples/Marc/material.config
View File

@ -1,429 +0,0 @@
#-------------------#
<homogenization>
#-------------------#
{../ConfigFiles/Homogenization_None_Dummy.config}
#-------------------#
<microstructure>
#-------------------#
[Grain001]
(constituent) phase 1 texture 1 fraction 1.0
[Grain002]
(constituent) phase 1 texture 2 fraction 1.0
[Grain003]
(constituent) phase 1 texture 3 fraction 1.0
[Grain004]
(constituent) phase 1 texture 4 fraction 1.0
[Grain005]
(constituent) phase 1 texture 5 fraction 1.0
[Grain006]
(constituent) phase 1 texture 6 fraction 1.0
[Grain007]
(constituent) phase 1 texture 7 fraction 1.0
[Grain008]
(constituent) phase 1 texture 8 fraction 1.0
[Grain009]
(constituent) phase 1 texture 9 fraction 1.0
[Grain010]
(constituent) phase 1 texture 10 fraction 1.0
[Grain011]
(constituent) phase 1 texture 11 fraction 1.0
[Grain012]
(constituent) phase 1 texture 12 fraction 1.0
[Grain013]
(constituent) phase 1 texture 13 fraction 1.0
[Grain014]
(constituent) phase 1 texture 14 fraction 1.0
[Grain015]
(constituent) phase 1 texture 15 fraction 1.0
[Grain016]
(constituent) phase 1 texture 16 fraction 1.0
[Grain017]
(constituent) phase 1 texture 17 fraction 1.0
[Grain018]
(constituent) phase 1 texture 18 fraction 1.0
[Grain019]
(constituent) phase 1 texture 19 fraction 1.0
[Grain020]
(constituent) phase 1 texture 20 fraction 1.0
[Grain021]
(constituent) phase 1 texture 21 fraction 1.0
[Grain022]
(constituent) phase 1 texture 22 fraction 1.0
[Grain023]
(constituent) phase 1 texture 23 fraction 1.0
[Grain024]
(constituent) phase 1 texture 24 fraction 1.0
[Grain025]
(constituent) phase 1 texture 25 fraction 1.0
[Grain026]
(constituent) phase 1 texture 26 fraction 1.0
[Grain027]
(constituent) phase 1 texture 27 fraction 1.0
[Grain028]
(constituent) phase 1 texture 28 fraction 1.0
[Grain029]
(constituent) phase 1 texture 29 fraction 1.0
[Grain030]
(constituent) phase 1 texture 30 fraction 1.0
[Grain031]
(constituent) phase 1 texture 31 fraction 1.0
[Grain032]
(constituent) phase 1 texture 32 fraction 1.0
[Grain033]
(constituent) phase 1 texture 33 fraction 1.0
[Grain034]
(constituent) phase 1 texture 34 fraction 1.0
[Grain035]
(constituent) phase 1 texture 35 fraction 1.0
[Grain036]
(constituent) phase 1 texture 36 fraction 1.0
[Grain037]
(constituent) phase 1 texture 37 fraction 1.0
[Grain038]
(constituent) phase 1 texture 38 fraction 1.0
[Grain039]
(constituent) phase 1 texture 39 fraction 1.0
[Grain040]
(constituent) phase 1 texture 40 fraction 1.0
[Grain041]
(constituent) phase 1 texture 41 fraction 1.0
[Grain042]
(constituent) phase 1 texture 42 fraction 1.0
[Grain043]
(constituent) phase 1 texture 43 fraction 1.0
[Grain044]
(constituent) phase 1 texture 44 fraction 1.0
[Grain045]
(constituent) phase 1 texture 45 fraction 1.0
[Grain046]
(constituent) phase 1 texture 46 fraction 1.0
[Grain047]
(constituent) phase 1 texture 47 fraction 1.0
[Grain048]
(constituent) phase 1 texture 48 fraction 1.0
[Grain049]
(constituent) phase 1 texture 49 fraction 1.0
[Grain050]
(constituent) phase 1 texture 50 fraction 1.0
[Grain051]
(constituent) phase 1 texture 51 fraction 1.0
[Grain052]
(constituent) phase 1 texture 52 fraction 1.0
[Grain053]
(constituent) phase 1 texture 53 fraction 1.0
[Grain054]
(constituent) phase 1 texture 54 fraction 1.0
[Grain055]
(constituent) phase 1 texture 55 fraction 1.0
[Grain056]
(constituent) phase 1 texture 56 fraction 1.0
[Grain057]
(constituent) phase 1 texture 57 fraction 1.0
[Grain058]
(constituent) phase 1 texture 58 fraction 1.0
[Grain059]
(constituent) phase 1 texture 59 fraction 1.0
[Grain060]
(constituent) phase 1 texture 60 fraction 1.0
[Grain061]
(constituent) phase 1 texture 61 fraction 1.0
[Grain062]
(constituent) phase 1 texture 62 fraction 1.0
[Grain063]
(constituent) phase 1 texture 63 fraction 1.0
[Grain064]
(constituent) phase 1 texture 64 fraction 1.0
[Grain065]
(constituent) phase 1 texture 65 fraction 1.0
[Grain066]
(constituent) phase 1 texture 66 fraction 1.0
[Grain067]
(constituent) phase 1 texture 67 fraction 1.0
[Grain068]
(constituent) phase 1 texture 68 fraction 1.0
[Grain069]
(constituent) phase 1 texture 69 fraction 1.0
[Grain070]
(constituent) phase 1 texture 70 fraction 1.0
[Grain071]
(constituent) phase 1 texture 71 fraction 1.0
[Grain072]
(constituent) phase 1 texture 72 fraction 1.0
[Grain073]
(constituent) phase 1 texture 73 fraction 1.0
[Grain074]
(constituent) phase 1 texture 74 fraction 1.0
[Grain075]
(constituent) phase 1 texture 75 fraction 1.0
[Grain076]
(constituent) phase 1 texture 76 fraction 1.0
[Grain077]
(constituent) phase 1 texture 77 fraction 1.0
[Grain078]
(constituent) phase 1 texture 78 fraction 1.0
[Grain079]
(constituent) phase 1 texture 79 fraction 1.0
[Grain080]
(constituent) phase 1 texture 80 fraction 1.0
[Grain081]
(constituent) phase 1 texture 81 fraction 1.0
[Grain082]
(constituent) phase 1 texture 82 fraction 1.0
[Grain083]
(constituent) phase 1 texture 83 fraction 1.0
[Grain084]
(constituent) phase 1 texture 84 fraction 1.0
[Grain085]
(constituent) phase 1 texture 85 fraction 1.0
[Grain086]
(constituent) phase 1 texture 86 fraction 1.0
[Grain087]
(constituent) phase 1 texture 87 fraction 1.0
[Grain088]
(constituent) phase 1 texture 88 fraction 1.0
[Grain089]
(constituent) phase 1 texture 89 fraction 1.0
[Grain090]
(constituent) phase 1 texture 90 fraction 1.0
[Grain091]
(constituent) phase 1 texture 91 fraction 1.0
[Grain092]
(constituent) phase 1 texture 92 fraction 1.0
[Grain093]
(constituent) phase 1 texture 93 fraction 1.0
[Grain094]
(constituent) phase 1 texture 94 fraction 1.0
[Grain095]
(constituent) phase 1 texture 95 fraction 1.0
[Grain096]
(constituent) phase 1 texture 96 fraction 1.0
[Grain097]
(constituent) phase 1 texture 97 fraction 1.0
[Grain098]
(constituent) phase 1 texture 98 fraction 1.0
[Grain099]
(constituent) phase 1 texture 99 fraction 1.0
[Grain100]
(constituent) phase 1 texture 100 fraction 1.0
[cubeGrain]
(constituent) phase 1 texture 101 fraction 1.0
#-------------------#
<texture>
#-------------------#
[Grain001]
(gauss) phi1 359.121452 Phi 82.319471 Phi2 347.729535
[Grain002]
(gauss) phi1 269.253967 Phi 105.379919 Phi2 173.029284
[Grain003]
(gauss) phi1 26.551535 Phi 171.606752 Phi2 124.949264
[Grain004]
(gauss) phi1 123.207774 Phi 124.339577 Phi2 47.937748
[Grain005]
(gauss) phi1 324.188825 Phi 103.089216 Phi2 160.373624
[Grain006]
(gauss) phi1 238.295585 Phi 165.416882 Phi2 234.307741
[Grain007]
(gauss) phi1 232.707177 Phi 110.733726 Phi2 308.049265
[Grain008]
(gauss) phi1 144.463291 Phi 125.891441 Phi2 348.674207
[Grain009]
(gauss) phi1 215.423832 Phi 69.759502 Phi2 164.477632
[Grain010]
(gauss) phi1 118.805444 Phi 143.057031 Phi2 271.963190
[Grain011]
(gauss) phi1 218.049576 Phi 64.017550 Phi2 323.040457
[Grain012]
(gauss) phi1 236.962483 Phi 134.312093 Phi2 220.433366
[Grain013]
(gauss) phi1 352.317686 Phi 3.356527 Phi2 92.447275
[Grain014]
(gauss) phi1 198.311545 Phi 71.452240 Phi2 199.441849
[Grain015]
(gauss) phi1 351.993635 Phi 36.500987 Phi2 236.852886
[Grain016]
(gauss) phi1 262.389063 Phi 101.249950 Phi2 334.305959
[Grain017]
(gauss) phi1 53.220668 Phi 69.570254 Phi2 277.061151
[Grain018]
(gauss) phi1 122.156119 Phi 140.207051 Phi2 221.172906
[Grain019]
(gauss) phi1 295.422170 Phi 26.595511 Phi2 263.206315
[Grain020]
(gauss) phi1 179.137406 Phi 104.500977 Phi2 151.742108
[Grain021]
(gauss) phi1 199.045094 Phi 5.228899 Phi2 356.542109
[Grain022]
(gauss) phi1 268.671476 Phi 24.835403 Phi2 33.578889
[Grain023]
(gauss) phi1 264.248527 Phi 59.766630 Phi2 340.865462
[Grain024]
(gauss) phi1 254.223491 Phi 51.125301 Phi2 201.094027
[Grain025]
(gauss) phi1 22.214008 Phi 92.248774 Phi2 215.168318
[Grain026]
(gauss) phi1 49.511491 Phi 79.933539 Phi2 187.188575
[Grain027]
(gauss) phi1 318.916204 Phi 113.102650 Phi2 241.076629
[Grain028]
(gauss) phi1 239.378433 Phi 89.578655 Phi2 94.167043
[Grain029]
(gauss) phi1 27.561421 Phi 142.892093 Phi2 197.735666
[Grain030]
(gauss) phi1 135.210581 Phi 165.859834 Phi2 285.449561
[Grain031]
(gauss) phi1 223.515916 Phi 56.824378 Phi2 343.289074
[Grain032]
(gauss) phi1 41.127974 Phi 111.289145 Phi2 214.855145
[Grain033]
(gauss) phi1 17.335045 Phi 140.496745 Phi2 77.747371
[Grain034]
(gauss) phi1 36.206421 Phi 148.574232 Phi2 88.870226
[Grain035]
(gauss) phi1 159.618336 Phi 125.680504 Phi2 204.119403
[Grain036]
(gauss) phi1 8.752464 Phi 99.173166 Phi2 143.227089
[Grain037]
(gauss) phi1 351.570753 Phi 67.343218 Phi2 1.779612
[Grain038]
(gauss) phi1 46.771572 Phi 155.018674 Phi2 302.319987
[Grain039]
(gauss) phi1 244.255976 Phi 80.566566 Phi2 264.069331
[Grain040]
(gauss) phi1 41.775388 Phi 47.109507 Phi2 300.598550
[Grain041]
(gauss) phi1 268.753103 Phi 46.654050 Phi2 190.382041
[Grain042]
(gauss) phi1 239.574480 Phi 62.517793 Phi2 147.817535
[Grain043]
(gauss) phi1 128.059775 Phi 61.916743 Phi2 169.674359
[Grain044]
(gauss) phi1 166.545156 Phi 58.709099 Phi2 252.885391
[Grain045]
(gauss) phi1 92.867691 Phi 28.906456 Phi2 164.197290
[Grain046]
(gauss) phi1 291.056147 Phi 35.145174 Phi2 250.155599
[Grain047]
(gauss) phi1 79.015862 Phi 44.772479 Phi2 267.982808
[Grain048]
(gauss) phi1 108.400702 Phi 69.883075 Phi2 222.737053
[Grain049]
(gauss) phi1 348.326500 Phi 11.339714 Phi2 121.682346
[Grain050]
(gauss) phi1 331.476209 Phi 108.775043 Phi2 335.139671
[Grain051]
(gauss) phi1 196.750278 Phi 93.955106 Phi2 63.689075
[Grain052]
(gauss) phi1 136.077875 Phi 130.508342 Phi2 128.468976
[Grain053]
(gauss) phi1 239.643513 Phi 76.284643 Phi2 168.821008
[Grain054]
(gauss) phi1 113.850670 Phi 117.531757 Phi2 71.971648
[Grain055]
(gauss) phi1 149.554071 Phi 16.543098 Phi2 195.556172
[Grain056]
(gauss) phi1 46.626579 Phi 52.447846 Phi2 304.495569
[Grain057]
(gauss) phi1 255.251821 Phi 86.678048 Phi2 238.982712
[Grain058]
(gauss) phi1 324.266133 Phi 28.075458 Phi2 41.191295
[Grain059]
(gauss) phi1 312.000332 Phi 74.648725 Phi2 87.403581
[Grain060]
(gauss) phi1 57.742481 Phi 163.241519 Phi2 68.491438
[Grain061]
(gauss) phi1 112.442447 Phi 51.735320 Phi2 206.538656
[Grain062]
(gauss) phi1 297.453842 Phi 115.283041 Phi2 57.785319
[Grain063]
(gauss) phi1 119.132681 Phi 117.923565 Phi2 196.121206
[Grain064]
(gauss) phi1 199.267314 Phi 163.091476 Phi2 53.549301
[Grain065]
(gauss) phi1 37.765215 Phi 76.795488 Phi2 146.264753
[Grain066]
(gauss) phi1 324.550183 Phi 27.665150 Phi2 56.383148
[Grain067]
(gauss) phi1 337.305377 Phi 136.807151 Phi2 133.661586
[Grain068]
(gauss) phi1 115.744041 Phi 64.536978 Phi2 262.694800
[Grain069]
(gauss) phi1 136.293403 Phi 48.862462 Phi2 343.319175
[Grain070]
(gauss) phi1 111.030931 Phi 80.823213 Phi2 84.041594
[Grain071]
(gauss) phi1 303.985249 Phi 118.929631 Phi2 302.307709
[Grain072]
(gauss) phi1 193.556259 Phi 75.928015 Phi2 176.696899
[Grain073]
(gauss) phi1 102.543259 Phi 121.929923 Phi2 234.496773
[Grain074]
(gauss) phi1 218.581323 Phi 101.753894 Phi2 305.566089
[Grain075]
(gauss) phi1 229.542114 Phi 118.839215 Phi2 129.179156
[Grain076]
(gauss) phi1 202.258840 Phi 139.205956 Phi2 352.248979
[Grain077]
(gauss) phi1 137.954289 Phi 63.806918 Phi2 128.975049
[Grain078]
(gauss) phi1 327.557366 Phi 84.987420 Phi2 345.483143
[Grain079]
(gauss) phi1 334.610243 Phi 74.535474 Phi2 106.419231
[Grain080]
(gauss) phi1 62.906243 Phi 46.752029 Phi2 222.692276
[Grain081]
(gauss) phi1 254.121439 Phi 121.005485 Phi2 287.265977
[Grain082]
(gauss) phi1 140.765045 Phi 141.268031 Phi2 271.327656
[Grain083]
(gauss) phi1 10.726984 Phi 66.339177 Phi2 189.073212
[Grain084]
(gauss) phi1 270.921536 Phi 72.821127 Phi2 313.590515
[Grain085]
(gauss) phi1 299.059668 Phi 23.884874 Phi2 80.016277
[Grain086]
(gauss) phi1 208.617406 Phi 11.031834 Phi2 302.388247
[Grain087]
(gauss) phi1 62.929967 Phi 65.223261 Phi2 108.558265
[Grain088]
(gauss) phi1 9.014959 Phi 33.542169 Phi2 247.970366
[Grain089]
(gauss) phi1 272.432808 Phi 30.065174 Phi2 19.803570
[Grain090]
(gauss) phi1 179.621980 Phi 151.763475 Phi2 61.871794
[Grain091]
(gauss) phi1 247.810321 Phi 112.752980 Phi2 264.668469
[Grain092]
(gauss) phi1 270.780630 Phi 102.037858 Phi2 31.602610
[Grain093]
(gauss) phi1 17.626672 Phi 56.032415 Phi2 245.079600
[Grain094]
(gauss) phi1 112.165186 Phi 87.390459 Phi2 182.086729
[Grain095]
(gauss) phi1 157.869381 Phi 79.905131 Phi2 107.037081
[Grain096]
(gauss) phi1 106.163846 Phi 148.477084 Phi2 350.980466
[Grain097]
(gauss) phi1 262.138550 Phi 58.923588 Phi2 111.303439
[Grain098]
(gauss) phi1 88.739397 Phi 119.092789 Phi2 222.502594
[Grain099]
(gauss) phi1 337.603765 Phi 10.145102 Phi2 80.934916
[Grain100]
(gauss) phi1 341.022242 Phi 45.927285 Phi2 252.045476
[cube]
(gauss) phi1 0 Phi 0 phi2 0
#-------------------#
<phase>
#-------------------#
{../ConfigFiles/Phase_Phenopowerlaw_Aluminum.config}
{../ConfigFiles/Phase_Isotropic_AluminumIsotropic.config}

31
examples/Marc/material.yaml Normal file
View File

@ -0,0 +1,31 @@
---
homogenization:
SX:
N_constituents: 1
mechanical: {type: pass}
phase:
Aluminum:
lattice: cF
mechanical:
output: [F, P, F_e, F_p, L_p, O]
elastic: {type: Hooke, C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9}
plastic:
type: phenopowerlaw
N_sl: [12]
a_sl: 2.25
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
xi_0_sl: [31e6]
xi_inf_sl: [63e6]
material:
- homogenization: SX
constituents:
- phase: Aluminum
v: 1.0
O: [0.9330127018922194, 0.25, 0.06698729810778066, 0.25]

BIN
examples/Marc/rotation.mud
View File

Binary file not shown.

406
examples/Marc/sheet_r-value.dat Normal file
View File

@ -0,0 +1,406 @@
title r-value
$....MARC input file produced by Marc Mentat 2019 (64bit)
$...................................
$....input file using extended precision
extended
$...................................
sizing 0 80 165 0
alloc 25
elements 7
version 14 1 0 1
table 0 0 2 1 1 0 0 1
processor 1 1 1 0
$no list
large stra 2 1 0 0 0 0 0
all points
no echo 1 2 3 4
state vars 3
end
$...................
solver
8 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
optimize 11
connectivity
0 0 1 0 1 1 0 0 0
1 7 2 5 20 17 1 4 19 16
2 7 3 6 21 18 2 5 20 17
3 7 5 8 23 20 4 7 22 19
4 7 6 9 24 21 5 8 23 20
5 7 8 11 26 23 7 10 25 22
6 7 9 12 27 24 8 11 26 23
7 7 11 14 29 26 10 13 28 25
8 7 12 15 30 27 11 14 29 26
9 7 17 20 35 32 16 19 34 31
10 7 18 21 36 33 17 20 35 32
11 7 20 23 38 35 19 22 37 34
12 7 21 24 39 36 20 23 38 35
13 7 23 26 41 38 22 25 40 37
14 7 24 27 42 39 23 26 41 38
15 7 26 29 44 41 25 28 43 40
16 7 27 30 45 42 26 29 44 41
17 7 32 35 50 47 31 34 49 46
18 7 33 36 51 48 32 35 50 47
19 7 35 38 53 50 34 37 52 49
20 7 36 39 54 51 35 38 53 50
21 7 38 41 56 53 37 40 55 52
22 7 39 42 57 54 38 41 56 53
23 7 41 44 59 56 40 43 58 55
24 7 42 45 60 57 41 44 59 56
25 7 47 50 65 62 46 49 64 61
26 7 48 51 66 63 47 50 65 62
27 7 50 53 68 65 49 52 67 64
28 7 51 54 69 66 50 53 68 65
29 7 53 56 71 68 52 55 70 67
30 7 54 57 72 69 53 56 71 68
31 7 56 59 74 71 55 58 73 70
32 7 57 60 75 72 56 59 74 71
33 7 62 65 80 77 61 64 79 76
34 7 63 66 81 78 62 65 80 77
35 7 65 68 83 80 64 67 82 79
36 7 66 69 84 81 65 68 83 80
37 7 68 71 86 83 67 70 85 82
38 7 69 72 87 84 68 71 86 83
39 7 71 74 89 86 70 73 88 85
40 7 72 75 90 87 71 74 89 86
41 7 77 80 95 92 76 79 94 91
42 7 78 81 96 93 77 80 95 92
43 7 80 83 98 95 79 82 97 94
44 7 81 84 99 96 80 83 98 95
45 7 83 86 101 98 82 85 100 97
46 7 84 87 102 99 83 86 101 98
47 7 86 89 104 101 85 88 103 100
48 7 87 90 105 102 86 89 104 101
49 7 92 95 110 107 91 94 109 106
50 7 93 96 111 108 92 95 110 107
51 7 95 98 113 110 94 97 112 109
52 7 96 99 114 111 95 98 113 110
53 7 98 101 116 113 97 100 115 112
54 7 99 102 117 114 98 101 116 113
55 7 101 104 119 116 100 103 118 115
56 7 102 105 120 117 101 104 119 116
57 7 107 110 125 122 106 109 124 121
58 7 108 111 126 123 107 110 125 122
59 7 110 113 128 125 109 112 127 124
60 7 111 114 129 126 110 113 128 125
61 7 113 116 131 128 112 115 130 127
62 7 114 117 132 129 113 116 131 128
63 7 116 119 134 131 115 118 133 130
64 7 117 120 135 132 116 119 134 131
65 7 122 125 140 137 121 124 139 136
66 7 123 126 141 138 122 125 140 137
67 7 125 128 143 140 124 127 142 139
68 7 126 129 144 141 125 128 143 140
69 7 128 131 146 143 127 130 145 142
70 7 129 132 147 144 128 131 146 143
71 7 131 134 149 146 130 133 148 145
72 7 132 135 150 147 131 134 149 146
73 7 137 140 155 152 136 139 154 151
74 7 138 141 156 153 137 140 155 152
75 7 140 143 158 155 139 142 157 154
76 7 141 144 159 156 140 143 158 155
77 7 143 146 161 158 142 145 160 157
78 7 144 147 162 159 143 146 161 158
79 7 146 149 164 161 145 148 163 160
80 7 147 150 165 162 146 149 164 161
coordinates
3 165 0 1
1-4.000000000000000+1-1.000000000000000+1-5.000000000000000-1
2-4.000000000000000+1-1.000000000000000+1 0.000000000000000+0
3-4.000000000000000+1-1.000000000000000+1 5.000000000000000-1
4-4.000000000000000+1-5.000000000000000+0-5.000000000000000-1
5-4.000000000000000+1-5.000000000000000+0 0.000000000000000+0
6-4.000000000000000+1-5.000000000000000+0 5.000000000000000-1
7-4.000000000000000+1 0.000000000000000+0-5.000000000000000-1
8-4.000000000000000+1 0.000000000000000+0 0.000000000000000+0
9-4.000000000000000+1 0.000000000000000+0 5.000000000000000-1
10-4.000000000000000+1 5.000000000000000+0-5.000000000000000-1
11-4.000000000000000+1 5.000000000000000+0 0.000000000000000+0
12-4.000000000000000+1 5.000000000000000+0 5.000000000000000-1
13-4.000000000000000+1 9.999999999999996+0-5.000000000000000-1
14-4.000000000000000+1 9.999999999999996+0 0.000000000000000+0
15-4.000000000000000+1 9.999999999999996+0 5.000000000000000-1
16-3.200000000000000+1-1.000000000000000+1-5.000000000000000-1
17-3.200000000000000+1-1.000000000000000+1 0.000000000000000+0
18-3.200000000000000+1-1.000000000000000+1 5.000000000000000-1
19-3.200000000000000+1-5.000000000000000+0-5.000000000000000-1
20-3.200000000000000+1-5.000000000000000+0 0.000000000000000+0
21-3.200000000000000+1-5.000000000000000+0 5.000000000000000-1
22-3.200000000000000+1 0.000000000000000+0-5.000000000000000-1
23-3.200000000000000+1 0.000000000000000+0 0.000000000000000+0
24-3.200000000000000+1 0.000000000000000+0 5.000000000000000-1
25-3.200000000000000+1 5.000000000000000+0-5.000000000000000-1
26-3.200000000000000+1 5.000000000000000+0 0.000000000000000+0
27-3.200000000000000+1 5.000000000000000+0 5.000000000000000-1
28-3.200000000000000+1 1.000000000000000+1-5.000000000000000-1
29-3.200000000000000+1 1.000000000000000+1 0.000000000000000+0
30-3.200000000000000+1 1.000000000000000+1 5.000000000000000-1
31-2.400000000000001+1-1.000000000000000+1-5.000000000000000-1
32-2.400000000000001+1-1.000000000000000+1 0.000000000000000+0
33-2.400000000000001+1-1.000000000000000+1 5.000000000000000-1
34-2.400000000000000+1-5.000000000000000+0-5.000000000000000-1
35-2.400000000000000+1-5.000000000000000+0 0.000000000000000+0
36-2.400000000000000+1-5.000000000000000+0 5.000000000000000-1
37-2.400000000000001+1 0.000000000000000+0-5.000000000000000-1
38-2.400000000000001+1 0.000000000000000+0 0.000000000000000+0
39-2.400000000000001+1 0.000000000000000+0 5.000000000000000-1
40-2.400000000000001+1 5.000000000000002+0-5.000000000000000-1
41-2.400000000000001+1 5.000000000000002+0 0.000000000000000+0
42-2.400000000000001+1 5.000000000000002+0 5.000000000000000-1
43-2.400000000000001+1 9.999999999999996+0-5.000000000000000-1
44-2.400000000000001+1 9.999999999999996+0 0.000000000000000+0
45-2.400000000000001+1 9.999999999999996+0 5.000000000000000-1
46-1.600000000000000+1-1.000000000000000+1-5.000000000000000-1
47-1.600000000000000+1-1.000000000000000+1 0.000000000000000+0
48-1.600000000000000+1-1.000000000000000+1 5.000000000000000-1
49-1.600000000000000+1-5.000000000000000+0-5.000000000000000-1
50-1.600000000000000+1-5.000000000000000+0 0.000000000000000+0
51-1.600000000000000+1-5.000000000000000+0 5.000000000000000-1
52-1.600000000000000+1 0.000000000000000+0-5.000000000000000-1
53-1.600000000000000+1 0.000000000000000+0 0.000000000000000+0
54-1.600000000000000+1 0.000000000000000+0 5.000000000000000-1
55-1.600000000000000+1 4.999999999999998+0-5.000000000000000-1
56-1.600000000000000+1 4.999999999999998+0 0.000000000000000+0
57-1.600000000000000+1 4.999999999999998+0 5.000000000000000-1
58-1.600000000000000+1 9.999999999999996+0-5.000000000000000-1
59-1.600000000000000+1 9.999999999999996+0 0.000000000000000+0
60-1.600000000000000+1 9.999999999999996+0 5.000000000000000-1
61-8.000000000000000+0-1.000000000000000+1-5.000000000000000-1
62-8.000000000000000+0-1.000000000000000+1 0.000000000000000+0
63-8.000000000000000+0-1.000000000000000+1 5.000000000000000-1
64-7.999999999999994+0-5.000000000000000+0-5.000000000000000-1
65-7.999999999999994+0-5.000000000000000+0 0.000000000000000+0
66-7.999999999999994+0-5.000000000000000+0 5.000000000000000-1
67-8.000000000000000+0 0.000000000000000+0-5.000000000000000-1
68-8.000000000000000+0 0.000000000000000+0 0.000000000000000+0
69-8.000000000000000+0 0.000000000000000+0 5.000000000000000-1
70-8.000000000000000+0 5.000000000000000+0-5.000000000000000-1
71-8.000000000000000+0 5.000000000000000+0 0.000000000000000+0
72-8.000000000000000+0 5.000000000000000+0 5.000000000000000-1
73-8.000000000000000+0 9.999999999999996+0-5.000000000000000-1
74-8.000000000000000+0 9.999999999999996+0 0.000000000000000+0
75-8.000000000000000+0 9.999999999999996+0 5.000000000000000-1
76 0.000000000000000+0-1.000000000000000+1-5.000000000000000-1
77 0.000000000000000+0-1.000000000000000+1 0.000000000000000+0
78 0.000000000000000+0-1.000000000000000+1 5.000000000000000-1
79 0.000000000000000+0-5.000000000000000+0-5.000000000000000-1
80 0.000000000000000+0-5.000000000000000+0 0.000000000000000+0
81 0.000000000000000+0-5.000000000000000+0 5.000000000000000-1
82 0.000000000000000+0 0.000000000000000+0-5.000000000000000-1
83 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0
84 0.000000000000000+0 0.000000000000000+0 5.000000000000000-1
85 0.000000000000000+0 5.000000000000000+0-5.000000000000000-1
86 0.000000000000000+0 5.000000000000000+0 0.000000000000000+0
87 0.000000000000000+0 5.000000000000000+0 5.000000000000000-1
88 0.000000000000000+0 9.999999999999996+0-5.000000000000000-1
89 0.000000000000000+0 9.999999999999996+0 0.000000000000000+0
90 0.000000000000000+0 9.999999999999996+0 5.000000000000000-1
91 8.000000000000000+0-1.000000000000000+1-5.000000000000000-1
92 8.000000000000000+0-1.000000000000000+1 0.000000000000000+0
93 8.000000000000000+0-1.000000000000000+1 5.000000000000000-1
94 8.000000000000000+0-5.000000000000000+0-5.000000000000000-1
95 8.000000000000000+0-5.000000000000000+0 0.000000000000000+0
96 8.000000000000000+0-5.000000000000000+0 5.000000000000000-1
97 8.000000000000000+0 0.000000000000000+0-5.000000000000000-1
98 8.000000000000000+0 0.000000000000000+0 0.000000000000000+0
99 8.000000000000000+0 0.000000000000000+0 5.000000000000000-1
100 8.000000000000000+0 5.000000000000000+0-5.000000000000000-1
101 8.000000000000000+0 5.000000000000000+0 0.000000000000000+0
102 8.000000000000000+0 5.000000000000000+0 5.000000000000000-1
103 7.999999999999994+0 9.999999999999996+0-5.000000000000000-1
104 7.999999999999994+0 9.999999999999996+0 0.000000000000000+0
105 7.999999999999994+0 9.999999999999996+0 5.000000000000000-1
106 1.600000000000000+1-1.000000000000000+1-5.000000000000000-1
107 1.600000000000000+1-1.000000000000000+1 0.000000000000000+0
108 1.600000000000000+1-1.000000000000000+1 5.000000000000000-1
109 1.600000000000000+1-5.000000000000000+0-5.000000000000000-1
110 1.600000000000000+1-5.000000000000000+0 0.000000000000000+0
111 1.600000000000000+1-5.000000000000000+0 5.000000000000000-1
112 1.600000000000000+1 0.000000000000000+0-5.000000000000000-1
113 1.600000000000000+1 0.000000000000000+0 0.000000000000000+0
114 1.600000000000000+1 0.000000000000000+0 5.000000000000000-1
115 1.600000000000000+1 5.000000000000000+0-5.000000000000000-1
116 1.600000000000000+1 5.000000000000000+0 0.000000000000000+0
117 1.600000000000000+1 5.000000000000000+0 5.000000000000000-1
118 1.600000000000000+1 9.999999999999996+0-5.000000000000000-1
119 1.600000000000000+1 9.999999999999996+0 0.000000000000000+0
120 1.600000000000000+1 9.999999999999996+0 5.000000000000000-1
121 2.400000000000000+1-1.000000000000000+1-5.000000000000000-1
122 2.400000000000000+1-1.000000000000000+1 0.000000000000000+0
123 2.400000000000000+1-1.000000000000000+1 5.000000000000000-1
124 2.400000000000001+1-5.000000000000000+0-5.000000000000000-1
125 2.400000000000001+1-5.000000000000000+0 0.000000000000000+0
126 2.400000000000001+1-5.000000000000000+0 5.000000000000000-1
127 2.400000000000000+1 0.000000000000000+0-5.000000000000000-1
128 2.400000000000000+1 0.000000000000000+0 0.000000000000000+0
129 2.400000000000000+1 0.000000000000000+0 5.000000000000000-1
130 2.400000000000001+1 5.000000000000000+0-5.000000000000000-1
131 2.400000000000001+1 5.000000000000000+0 0.000000000000000+0
132 2.400000000000001+1 5.000000000000000+0 5.000000000000000-1
133 2.400000000000000+1 9.999999999999996+0-5.000000000000000-1
134 2.400000000000000+1 9.999999999999996+0 0.000000000000000+0
135 2.400000000000000+1 9.999999999999996+0 5.000000000000000-1
136 3.200000000000000+1-1.000000000000000+1-5.000000000000000-1
137 3.200000000000000+1-1.000000000000000+1 0.000000000000000+0
138 3.200000000000000+1-1.000000000000000+1 5.000000000000000-1
139 3.199999999999999+1-5.000000000000000+0-5.000000000000000-1
140 3.199999999999999+1-5.000000000000000+0 0.000000000000000+0
141 3.199999999999999+1-5.000000000000000+0 5.000000000000000-1
142 3.200000000000000+1 0.000000000000000+0-5.000000000000000-1
143 3.200000000000000+1 0.000000000000000+0 0.000000000000000+0
144 3.200000000000000+1 0.000000000000000+0 5.000000000000000-1
145 3.200000000000000+1 5.000000000000000+0-5.000000000000000-1
146 3.200000000000000+1 5.000000000000000+0 0.000000000000000+0
147 3.200000000000000+1 5.000000000000000+0 5.000000000000000-1
148 3.199999999999999+1 9.999999999999996+0-5.000000000000000-1
149 3.199999999999999+1 9.999999999999996+0 0.000000000000000+0
150 3.199999999999999+1 9.999999999999996+0 5.000000000000000-1
151 3.999999999999999+1-1.000000000000000+1-5.000000000000000-1
152 3.999999999999999+1-1.000000000000000+1 0.000000000000000+0
153 3.999999999999999+1-1.000000000000000+1 5.000000000000000-1
154 3.999999999999999+1-5.000000000000000+0-5.000000000000000-1
155 3.999999999999999+1-5.000000000000000+0 0.000000000000000+0
156 3.999999999999999+1-5.000000000000000+0 5.000000000000000-1
157 3.999999999999999+1 0.000000000000000+0-5.000000000000000-1
158 3.999999999999999+1 0.000000000000000+0 0.000000000000000+0
159 3.999999999999999+1 0.000000000000000+0 5.000000000000000-1
160 3.999999999999999+1 5.000000000000000+0-5.000000000000000-1
161 3.999999999999999+1 5.000000000000000+0 0.000000000000000+0
162 3.999999999999999+1 5.000000000000000+0 5.000000000000000-1
163 3.999999999999999+1 9.999999999999996+0-5.000000000000000-1
164 3.999999999999999+1 9.999999999999996+0 0.000000000000000+0
165 3.999999999999999+1 9.999999999999996+0 5.000000000000000-1
define element set Material_Nummer_elements
1 to 80
define node set unten_y_nodes
2 5 8 11 14
define node set oben_y_nodes
152 155 158 161 164
define node set unten_fest_nodes
1 to 15
define node set oben_ziehen_nodes
151 to 165
define node set unten_z_nodes
7 to 9
define node set oben_z_nodes
157 to 159
define element set texture_elements
1 to 80
hypoelastic
1 0 1 0 1TKS 0
1.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0
0 0 0 0 0 0 0
mat color
1 1 230 0 0
table weg_x
1 1 0 0 2
1 2 2 0 0 2 0 0 2 0 0 2
0.000000000000000+0 0.000000000000000+0
2.000000000000000+2 1.600000000000000+1
geometry
0 0 2
1 9 1 230 0 0
r-value-sample
0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0
usdata 1
fixed disp
1 0 0 0 1 0unten_z
0.000000000000000+0 0.000000000000000+0
0 0
1 2
2
unten_z_nodes
1 0 0 0 1 0unten_y
0.000000000000000+0 0.000000000000000+0
0 0
1 3
2
unten_y_nodes
1 0 0 0 1 0oben_z
1.000000000000000+0 0.000000000000000+0
1 0
1 2
2
oben_z_nodes
1 0 0 0 1 0oben_y
1.000000000000000+0 0.000000000000000+0
1 0
1 3
2
oben_y_nodes
1 0 0 0 1 0unten_fest
0.000000000000000+0
0
1
2
unten_fest_nodes
1 0 0 0 1 0oben_ziehen
1.000000000000000+0
1
1
2
oben_ziehen_nodes
initial state
2 6 1 0 0 0Material_Nummer
1.000000000000000+0
0
1
Material_Nummer_elements
initial state
3 6 1 0 0 0texture
1.000000000000000+0
0
1
texture_elements
loadcase r-value
5
Material_Nummer
texture
unten_z
unten_y
unten_fest
no print
post
6 16 17 0 0 19 20 0 1 0 0 0 0 0 0 0
parameters
1.000000000000000+0 1.000000000000000+9 1.000000000000000+2 1.000000000000000+6 2.500000000000000-1 5.000000000000000-1 1.500000000000000+0-5.000000000000000-1
8.625000000000000+0 2.000000000000000+1 1.000000000000000-4 1.000000000000000-6 1.000000000000000+0 1.000000000000000-4
8.314000000000000+0 2.731500000000000+2 5.000000000000000-1 0.000000000000000+0 5.670510000000000-8 1.438769000000000-2 2.997900000000000+8 1.00000000000000+30
0.000000000000000+0 0.000000000000000+0 1.000000000000000+2 0.000000000000000+0 1.000000000000000+0-2.000000000000000+0 1.000000000000000+6 3.000000000000000+0
0.000000000000000+0 0.000000000000000+0 1.256637061000000-6 8.85418781700000-12 1.200000000000000+2 1.000000000000000-3 1.600000000000000+2 0.000000000000000+0
3.000000000000000+0 4.000000000000000-1
end option
$...................
$....start of loadcase Tensile
title Tensile
loadcase Tensile
6
unten_z
unten_y
oben_z
oben_y
unten_fest
oben_ziehen
control
99999 10 0 0 0 1 0 0 1 0 0 0 0 0 0
1.000000000000000-1 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0-1.000000000000000+0 0.000000000000000+0
parameters
1.000000000000000+0 1.000000000000000+9 1.000000000000000+2 1.000000000000000+6 2.500000000000000-1 5.000000000000000-1 1.500000000000000+0-5.000000000000000-1
8.625000000000000+0 2.000000000000000+1 1.000000000000000-4 1.000000000000000-6 1.000000000000000+0 1.000000000000000-4
8.314000000000000+0 2.731500000000000+2 5.000000000000000-1 0.000000000000000+0 5.670510000000000-8 1.438769000000000-2 2.997900000000000+8 1.00000000000000+30
0.000000000000000+0 0.000000000000000+0 1.000000000000000+2 0.000000000000000+0 1.000000000000000+0-1.000000000000000+0 1.000000000000000+6 3.000000000000000+0
0.000000000000000+0 0.000000000000000+0 1.256637061000000-6 8.85418781700000-12 1.200000000000000+2 1.000000000000000-3 1.600000000000000+2 0.000000000000000+0
3.000000000000000+0 4.000000000000000-1
auto load
100 0 10 0 0
time step
2.000000000000000+0
continue
$....end of loadcase Tensile
$...................

31
python/damask/solver/_marc.py
View File

@ -1,7 +1,6 @@
import subprocess import subprocess
import shlex import shlex
import re import re
import io
import os import os
from pathlib import Path from pathlib import Path
@ -41,13 +40,9 @@ class Marc:
return path_tools return path_tools
def submit_job(self, def submit_job(self, model, job,
model,
job = 'job1',
logfile = False,
compile = False, compile = False,
optimization = '', optimization = ''):
):
usersub = Path(os.environ['DAMASK_ROOT'])/'src/DAMASK_Marc' usersub = Path(os.environ['DAMASK_ROOT'])/'src/DAMASK_Marc'
usersub = usersub.parent/(usersub.name + ('.f90' if compile else '.marc')) usersub = usersub.parent/(usersub.name + ('.f90' if compile else '.marc'))
@ -64,22 +59,16 @@ class Marc:
' -prog ' + str(usersub.with_suffix('')) ' -prog ' + str(usersub.with_suffix(''))
print(cmd) print(cmd)
if logfile is not None: ret = subprocess.run(shlex.split(cmd),capture_output=True)
try:
f = open(logfile,'w+',newline='\n')
except TypeError:
f = logfile
else:
f = io.StringIO()
proc = subprocess.Popen(shlex.split(cmd),stdout=f,stderr=subprocess.STDOUT)
proc.wait()
f.seek(0)
try: try:
v = int(re.search('Exit number ([0-9]+)',''.join(f.readlines())).group(1)) v = int(re.search('Exit number ([0-9]+)',ret.stderr.decode()).group(1))
if 3004 != v:
print(ret.stderr.decode())
print(ret.stdout.decode())
raise RuntimeError(f'Marc simulation failed ({v})')
except (AttributeError,ValueError): except (AttributeError,ValueError):
print(ret.stderr.decode())
print(ret.stdout.decode())
raise RuntimeError('Marc simulation failed (unknown return value)') raise RuntimeError('Marc simulation failed (unknown return value)')
if v != 3004:
raise RuntimeError(f'Marc simulation failed ({v})')

21
src/CPFEM.f90
View File

@ -150,15 +150,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
H H
integer(pInt) elCP, & ! crystal plasticity element number integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog, mySource,ma i, j, k, l, m, n, ph, homog, mySource,ce
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll
elCP = mesh_FEM2DAMASK_elem(elFE) elCP = discretization_Marc_FEM2DAMASK_elem(elFE)
ce = discretization_Marc_FEM2DAMASK_cell(ip,elFE)
ma = (elCP-1) * discretization_nIPs + ip
if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then
print'(/,a)', '#############################################' print'(/,a)', '#############################################'
@ -183,8 +182,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
! temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & ! temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
! temperature_inp ! temperature_inp
!end select chosenThermal1 !end select chosenThermal1
homogenization_F0(1:3,1:3,ma) = ffn homogenization_F0(1:3,1:3,ce) = ffn
homogenization_F(1:3,1:3,ma) = ffn1 homogenization_F(1:3,1:3,ce) = ffn1
if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then
@ -209,17 +208,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else terminalIllness else terminalIllness
! translate from P to sigma ! translate from P to sigma
Kirchhoff = matmul(homogenization_P(1:3,1:3,ma), transpose(homogenization_F(1:3,1:3,ma))) Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(homogenization_F(1:3,1:3,ce)))
J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ma)) J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ce))
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
! translate from dP/dF to dCS/dE ! translate from dP/dF to dCS/dE
H = 0.0_pReal H = 0.0_pReal
do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
H(i,j,k,l) = H(i,j,k,l) & H(i,j,k,l) = H(i,j,k,l) &
+ homogenization_F(j,m,ma) * homogenization_F(l,n,ma) & + homogenization_F(j,m,ce) * homogenization_F(l,n,ce) &
* homogenization_dPdF(i,m,k,n,ma) & * homogenization_dPdF(i,m,k,n,ce) &
- math_delta(j,l) * homogenization_F(i,m,ma) * homogenization_P(k,m,ma) & - math_delta(j,l) * homogenization_F(i,m,ce) * homogenization_P(k,m,ce) &
+ 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) & + 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
+ Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k)) + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo

15
src/DAMASK_Marc.f90
View File

@ -218,7 +218,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
logical :: cutBack logical :: cutBack
real(pReal), dimension(6) :: stress real(pReal), dimension(6) :: stress
real(pReal), dimension(6,6) :: ddsdde real(pReal), dimension(6,6) :: ddsdde
integer :: computationMode, i, cp_en, node, CPnodeID integer :: computationMode, i, node, CPnodeID
integer(pI32) :: defaultNumThreadsInt !< default value set by Marc integer(pI32) :: defaultNumThreadsInt !< default value set by Marc
integer(pInt), save :: & integer(pInt), save :: &
@ -266,7 +266,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
computationMode = CPFEM_RESTOREJACOBIAN computationMode = CPFEM_RESTOREJACOBIAN
elseif (lovl == 6) then ! stress requested by marc elseif (lovl == 6) then ! stress requested by marc
computationMode = CPFEM_CALCRESULTS computationMode = CPFEM_CALCRESULTS
cp_en = mesh_FEM2DAMASK_elem(m(1))
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
terminallyIll = .false. terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0 cycleCounter = -1 ! first calc step increments this to cycle = 0
@ -344,8 +343,8 @@ subroutine flux(f,ts,n,time)
real(pReal), dimension(2), intent(out) :: & real(pReal), dimension(2), intent(out) :: &
f f
f(1) = homogenization_f_T(discretization_Marc_FEM2DAMASK_cell(n(3),n(1)))
f(2) = 0.0_pReal f(2) = 0.0_pReal
call thermal_conduction_getSource(f(1), n(3),mesh_FEM2DAMASK_elem(n(1)))
end subroutine flux end subroutine flux
@ -370,11 +369,11 @@ subroutine uedinc(inc,incsub)
if (inc > inc_written) then if (inc > inc_written) then
allocate(d_n(3,count(mesh_FEM2DAMASK_node /= -1))) allocate(d_n(3,count(discretization_Marc_FEM2DAMASK_node /= -1)))
do n = lbound(mesh_FEM2DAMASK_node,1), ubound(mesh_FEM2DAMASK_node,1) do n = lbound(discretization_Marc_FEM2DAMASK_node,1), ubound(discretization_Marc_FEM2DAMASK_node,1)
if (mesh_FEM2DAMASK_node(n) /= -1) then if (discretization_Marc_FEM2DAMASK_node(n) /= -1) then
call nodvar(1,n,d_n(1:3,mesh_FEM2DAMASK_node(n)),nqncomp,nqdatatype) call nodvar(1,n,d_n(1:3,discretization_Marc_FEM2DAMASK_node(n)),nqncomp,nqdatatype)
if(nqncomp == 2) d_n(3,mesh_FEM2DAMASK_node(n)) = 0.0_pReal if(nqncomp == 2) d_n(3,discretization_Marc_FEM2DAMASK_node(n)) = 0.0_pReal
endif endif
enddo enddo

56
src/Marc/discretization_Marc.f90
View File

@ -24,8 +24,8 @@ module discretization_marc
mesh_unitlength !< physical length of one unit in mesh MD: needs systematic_name mesh_unitlength !< physical length of one unit in mesh MD: needs systematic_name
integer, dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID MD: Needs systematic name discretization_Marc_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID
mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID MD: needs systematic_name discretization_Marc_FEM2DAMASK_node !< DAMASK node ID for Marc node ID
type tCellNodeDefinition type tCellNodeDefinition
@ -39,8 +39,9 @@ module discretization_marc
connectivity_cell !< cell connectivity for each element,ip/cell connectivity_cell !< cell connectivity for each element,ip/cell
public :: & public :: &
discretization_marc_init, & discretization_Marc_init, &
discretization_marc_UpdateNodeAndIpCoords discretization_Marc_updateNodeAndIpCoords, &
discretization_Marc_FEM2DAMASK_cell
contains contains
@ -48,7 +49,7 @@ contains
!> @brief initializes the mesh by calling all necessary private routines the mesh module !> @brief initializes the mesh by calling all necessary private routines the mesh module
!! Order and routines strongly depend on type of solver !! Order and routines strongly depend on type of solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine discretization_marc_init subroutine discretization_Marc_init
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
node0_elem, & !< node x,y,z coordinates (initially!) node0_elem, & !< node x,y,z coordinates (initially!)
@ -96,7 +97,7 @@ subroutine discretization_marc_init
call buildCells(connectivity_cell,cellNodeDefinition,& call buildCells(connectivity_cell,cellNodeDefinition,&
elem,connectivity_elem) elem,connectivity_elem)
node0_cell = buildCellNodes(node0_elem) node0_cell = buildCellNodes(node0_elem)
IP_reshaped = buildIPcoordinates(node0_cell) IP_reshaped = buildIPcoordinates(node0_cell)
call discretization_init(materialAt, IP_reshaped, node0_cell) call discretization_init(materialAt, IP_reshaped, node0_cell)
@ -114,25 +115,42 @@ subroutine discretization_marc_init
call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(elem)) call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(elem))
call geometry_plastic_nonlocal_results call geometry_plastic_nonlocal_results
end subroutine discretization_marc_init end subroutine discretization_Marc_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate and set current nodal and IP positions (including cell nodes) !> @brief Calculate and set current nodal and IP positions (including cell nodes)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine discretization_marc_UpdateNodeAndIpCoords(d_n) subroutine discretization_Marc_updateNodeAndIpCoords(d_n)
real(pReal), dimension(:,:), intent(in) :: d_n real(pReal), dimension(:,:), intent(in) :: d_n
real(pReal), dimension(:,:), allocatable :: node_cell real(pReal), dimension(:,:), allocatable :: node_cell
node_cell = buildCellNodes(discretization_NodeCoords0(1:3,1:maxval(mesh_FEM2DAMASK_node)) + d_n) node_cell = buildCellNodes(discretization_NodeCoords0(1:3,1:maxval(discretization_Marc_FEM2DAMASK_node)) + d_n)
call discretization_setNodeCoords(node_cell) call discretization_setNodeCoords(node_cell)
call discretization_setIPcoords(buildIPcoordinates(node_cell)) call discretization_setIPcoords(buildIPcoordinates(node_cell))
end subroutine discretization_marc_UpdateNodeAndIpCoords end subroutine discretization_Marc_updateNodeAndIpCoords
!--------------------------------------------------------------------------------------------------
!> @brief Calculate and set current nodal and IP positions (including cell nodes)
!--------------------------------------------------------------------------------------------------
function discretization_marc_FEM2DAMASK_cell(IP_FEM,elem_FEM) result(cell)
integer, intent(in) :: IP_FEM, elem_FEM
integer :: cell
real(pReal), dimension(:,:), allocatable :: node_cell
cell = (discretization_Marc_FEM2DAMASK_elem(elem_FEM)-1)*discretization_nIPs + IP_FEM
end function discretization_marc_FEM2DAMASK_cell
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -219,10 +237,10 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,materialAt)
call inputRead_elemType(elem, & call inputRead_elemType(elem, &
nElems,inputFile) nElems,inputFile)
call inputRead_mapElems(mesh_FEM2DAMASK_elem,& call inputRead_mapElems(discretization_Marc_FEM2DAMASK_elem,&
nElems,elem%nNodes,inputFile) nElems,elem%nNodes,inputFile)
call inputRead_mapNodes(mesh_FEM2DAMASK_node,& call inputRead_mapNodes(discretization_Marc_FEM2DAMASK_node,&
nNodes,inputFile) nNodes,inputFile)
call inputRead_elemNodes(node0_elem, & call inputRead_elemNodes(node0_elem, &
@ -532,7 +550,7 @@ subroutine inputRead_elemNodes(nodes, &
if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then
chunkPos = [4,1,10,11,30,31,50,51,70] chunkPos = [4,1,10,11,30,31,50,51,70]
do i=1,nNode do i=1,nNode
m = mesh_FEM2DAMASK_node(IO_intValue(fileContent(l+1+i),chunkPos,1)) m = discretization_Marc_FEM2DAMASK_node(IO_intValue(fileContent(l+1+i),chunkPos,1))
do j = 1,3 do j = 1,3
nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1) nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1)
enddo enddo
@ -653,11 +671,11 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
j = 0 j = 0
do i = 1,nElem do i = 1,nElem
chunkPos = IO_stringPos(fileContent(l+1+i+j)) chunkPos = IO_stringPos(fileContent(l+1+i+j))
e = mesh_FEM2DAMASK_elem(IO_intValue(fileContent(l+1+i+j),chunkPos,1)) e = discretization_Marc_FEM2DAMASK_elem(IO_intValue(fileContent(l+1+i+j),chunkPos,1))
if (e /= 0) then ! disregard non CP elems if (e /= 0) then ! disregard non CP elems
do k = 1,chunkPos(1)-2 do k = 1,chunkPos(1)-2
inputRead_connectivityElem(k,e) = & inputRead_connectivityElem(k,e) = &
mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2)) discretization_Marc_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2))
enddo enddo
nNodesAlreadyRead = chunkPos(1) - 2 nNodesAlreadyRead = chunkPos(1) - 2
do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line
@ -665,7 +683,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
chunkPos = IO_stringPos(fileContent(l+1+i+j)) chunkPos = IO_stringPos(fileContent(l+1+i+j))
do k = 1,chunkPos(1) do k = 1,chunkPos(1)
inputRead_connectivityElem(nNodesAlreadyRead+k,e) = & inputRead_connectivityElem(nNodesAlreadyRead+k,e) = &
mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k)) discretization_Marc_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k))
enddo enddo
nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
enddo enddo
@ -718,7 +736,7 @@ subroutine inputRead_material(materialAt,&
if (initialcondTableStyle == 2) m = m + 2 if (initialcondTableStyle == 2) m = m + 2
contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements
do i = 1,contInts(1) do i = 1,contInts(1)
e = mesh_FEM2DAMASK_elem(contInts(1+i)) e = discretization_Marc_FEM2DAMASK_elem(contInts(1+i))
materialAt(e) = myVal materialAt(e) = myVal
enddo enddo
if (initialcondTableStyle == 0) m = m + 1 if (initialcondTableStyle == 0) m = m + 1

2
src/grid/discretization_grid.f90
View File

@ -74,8 +74,10 @@ subroutine discretization_grid_init(restart)
allocate(materialAt_global(0)) ! needed for IntelMPI allocate(materialAt_global(0)) ! needed for IntelMPI
endif endif
call MPI_Bcast(grid,3,MPI_INTEGER,0,PETSC_COMM_WORLD, ierr) call MPI_Bcast(grid,3,MPI_INTEGER,0,PETSC_COMM_WORLD, ierr)
if (ierr /= 0) error stop 'MPI error' if (ierr /= 0) error stop 'MPI error'
if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1')
call MPI_Bcast(geomSize,3,MPI_DOUBLE,0,PETSC_COMM_WORLD, ierr) call MPI_Bcast(geomSize,3,MPI_DOUBLE,0,PETSC_COMM_WORLD, ierr)
if (ierr /= 0) error stop 'MPI error' if (ierr /= 0) error stop 'MPI error'
call MPI_Bcast(origin,3,MPI_DOUBLE,0,PETSC_COMM_WORLD, ierr) call MPI_Bcast(origin,3,MPI_DOUBLE,0,PETSC_COMM_WORLD, ierr)

48
src/grid/grid_damage_spectral.f90
View File

@ -15,8 +15,8 @@ module grid_damage_spectral
use IO use IO
use spectral_utilities use spectral_utilities
use discretization_grid use discretization_grid
use YAML_types
use homogenization use homogenization
use YAML_types
use config use config
implicit none implicit none
@ -61,7 +61,7 @@ contains
!> @brief allocates all neccessary fields and fills them with data !> @brief allocates all neccessary fields and fills them with data
! ToDo: Restart not implemented ! ToDo: Restart not implemented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_init subroutine grid_damage_spectral_init()
PetscInt, dimension(0:worldsize-1) :: localK PetscInt, dimension(0:worldsize-1) :: localK
DM :: damage_grid DM :: damage_grid
@ -88,10 +88,10 @@ subroutine grid_damage_spectral_init
num_generic => config_numerics%get('generic',defaultVal=emptyDict) num_generic => config_numerics%get('generic',defaultVal=emptyDict)
num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal)
if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
if (num%eps_damage_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_atol') if (num%eps_damage_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_atol')
if (num%eps_damage_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_rtol') if (num%eps_damage_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_rtol')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
@ -146,6 +146,7 @@ subroutine grid_damage_spectral_init
allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr) call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr)
call updateReference call updateReference
@ -216,21 +217,21 @@ end function grid_damage_spectral_solution
subroutine grid_damage_spectral_forward(cutBack) subroutine grid_damage_spectral_forward(cutBack)
logical, intent(in) :: cutBack logical, intent(in) :: cutBack
integer :: i, j, k, ce integer :: i, j, k, ce
DM :: dm_local DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr PetscErrorCode :: ierr
if (cutBack) then if (cutBack) then
phi_current = phi_lastInc phi_current = phi_lastInc
phi_stagInc = phi_lastInc phi_stagInc = phi_lastInc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reverting damage field state ! reverting damage field state
ce = 0
call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call homogenization_set_phi(phi_current(i,j,k),ce) call homogenization_set_phi(phi_current(i,j,k),ce)
@ -259,7 +260,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: i, j, k, ce integer :: i, j, k, ce
real(pReal) :: phiDot, mobility
phi_current = x_scal phi_current = x_scal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -272,8 +272,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(ce) - K_ref, & vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k))
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
@ -281,10 +280,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call damage_nonlocal_getSourceAndItsTangent(phiDot, phi_current(i,j,k),ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
mobility = damage_nonlocal_getMobility(ce) + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) &
+ mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
+ mu_ref*phi_current(i,j,k) + mu_ref*phi_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -293,6 +290,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
call utilities_FFTscalarForward call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc) call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc)
call utilities_FFTscalarBackward call utilities_FFTscalarBackward
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) & where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) & where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) &
@ -308,18 +306,18 @@ end subroutine formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief update reference viscosity and conductivity !> @brief update reference viscosity and conductivity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine updateReference subroutine updateReference()
integer :: ce,ierr
integer :: i,j,k,ce,ierr
ce = 0
K_ref = 0.0_pReal K_ref = 0.0_pReal
mu_ref = 0.0_pReal mu_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do ce = 1, product(grid(1:2))*grid3
ce = ce + 1 K_ref = K_ref + homogenization_K_phi(ce)
K_ref = K_ref + damage_nonlocal_getDiffusion(ce) mu_ref = mu_ref + homogenization_mu_phi(ce)
mu_ref = mu_ref + damage_nonlocal_getMobility(ce) enddo
enddo; enddo; enddo
K_ref = K_ref*wgt K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mu_ref = mu_ref*wgt mu_ref = mu_ref*wgt

54
src/grid/grid_thermal_spectral.f90
View File

@ -18,7 +18,6 @@ module grid_thermal_spectral
use homogenization use homogenization
use YAML_types use YAML_types
use config use config
use material
implicit none implicit none
private private
@ -46,7 +45,7 @@ module grid_thermal_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc. ! reference diffusion tensor, mobility etc.
integer :: totalIter = 0 !< total iteration in current increment integer :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3) :: K_ref real(pReal), dimension(3,3) :: K_ref
real(pReal) :: mu_ref real(pReal) :: mu_ref
@ -68,7 +67,7 @@ subroutine grid_thermal_spectral_init(T_0)
PetscInt, dimension(0:worldsize-1) :: localK PetscInt, dimension(0:worldsize-1) :: localK
integer :: i, j, k, ce integer :: i, j, k, ce
DM :: thermal_grid DM :: thermal_grid
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: T_PETSc
PetscErrorCode :: ierr PetscErrorCode :: ierr
class(tNode), pointer :: & class(tNode), pointer :: &
num_grid num_grid
@ -85,9 +84,9 @@ subroutine grid_thermal_spectral_init(T_0)
num%eps_thermal_atol = num_grid%get_asFloat ('eps_thermal_atol',defaultVal=1.0e-2_pReal) num%eps_thermal_atol = num_grid%get_asFloat ('eps_thermal_atol',defaultVal=1.0e-2_pReal)
num%eps_thermal_rtol = num_grid%get_asFloat ('eps_thermal_rtol',defaultVal=1.0e-6_pReal) num%eps_thermal_rtol = num_grid%get_asFloat ('eps_thermal_rtol',defaultVal=1.0e-6_pReal)
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
if (num%eps_thermal_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_thermal_atol') if (num%eps_thermal_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_thermal_atol')
if (num%eps_thermal_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_thermal_rtol') if (num%eps_thermal_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_thermal_rtol')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
@ -130,6 +129,7 @@ subroutine grid_thermal_spectral_init(T_0)
allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal) allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
@ -138,9 +138,10 @@ subroutine grid_thermal_spectral_init(T_0)
T_stagInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k)
call homogenization_thermal_setField(T_0,0.0_pReal,ce) call homogenization_thermal_setField(T_0,0.0_pReal,ce)
enddo; enddo; enddo enddo; enddo; enddo
call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) T_PETSc(xstart:xend,ystart:yend,zstart:zend) = T_current
call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,T_PETSc,ierr); CHKERRQ(ierr)
call updateReference call updateReference
@ -190,9 +191,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k), & call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce)
(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, &
ce)
enddo; enddo; enddo enddo; enddo; enddo
call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
@ -223,16 +222,14 @@ subroutine grid_thermal_spectral_forward(cutBack)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reverting thermal field state ! reverting thermal field state
ce = 0
call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr) call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k), & call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce)
(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, &
ce)
enddo; enddo; enddo enddo; enddo; enddo
else else
T_lastInc = T_current T_lastInc = T_current
@ -258,7 +255,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: i, j, k, ce integer :: i, j, k, ce
real(pReal) :: Tdot
T_current = x_scal T_current = x_scal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -271,8 +267,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(ce) - K_ref, & vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k))
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field
@ -280,9 +275,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call thermal_conduction_getSource(Tdot,1,ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_T(ce)) &
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
+ homogenization_thermal_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
+ mu_ref*T_current(i,j,k) + mu_ref*T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -302,18 +296,18 @@ end subroutine formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief update reference viscosity and conductivity !> @brief update reference viscosity and conductivity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine updateReference subroutine updateReference()
integer :: ce,ierr
integer :: i,j,k,ce,ierr
ce = 0
K_ref = 0.0_pReal K_ref = 0.0_pReal
mu_ref = 0.0_pReal mu_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do ce = 1, product(grid(1:2))*grid3
ce = ce + 1 K_ref = K_ref + homogenization_K_T(ce)
K_ref = K_ref + thermal_conduction_getConductivity(ce) mu_ref = mu_ref + homogenization_mu_T(ce)
mu_ref = mu_ref + homogenization_thermal_mu_T(ce) enddo
enddo; enddo; enddo
K_ref = K_ref*wgt K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mu_ref = mu_ref*wgt mu_ref = mu_ref*wgt

91
src/homogenization.f90
View File

@ -100,10 +100,6 @@ module homogenization
integer, intent(in) :: ce integer, intent(in) :: ce
end subroutine damage_partition end subroutine damage_partition
module subroutine thermal_homogenize(ip,el)
integer, intent(in) :: ip,el
end subroutine thermal_homogenize
module subroutine mechanical_homogenize(dt,ce) module subroutine mechanical_homogenize(dt,ce)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
@ -135,46 +131,41 @@ module homogenization
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
end function mechanical_updateState end function mechanical_updateState
module function homogenization_mu_T(ce) result(mu)
integer, intent(in) :: ce
real(pReal) :: mu
end function homogenization_mu_T
module function thermal_conduction_getConductivity(ce) result(K) module function homogenization_K_T(ce) result(K)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), dimension(3,3) :: K real(pReal), dimension(3,3) :: K
end function thermal_conduction_getConductivity end function homogenization_K_T
module function homogenization_thermal_mu_T(ce) result(mu_T) module function homogenization_f_T(ce) result(f)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: mu_T real(pReal) :: f
end function homogenization_thermal_mu_T end function homogenization_f_T
module subroutine homogenization_thermal_setField(T,dot_T, ce) module subroutine homogenization_thermal_setField(T,dot_T, ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
end subroutine homogenization_thermal_setField end subroutine homogenization_thermal_setField
module function homogenization_thermal_T(ce) result(T) module function homogenization_mu_phi(ce) result(mu)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: T real(pReal) :: mu
end function homogenization_thermal_T end function homogenization_mu_phi
module subroutine thermal_conduction_getSource(Tdot, ip, el) module function homogenization_K_phi(ce) result(K)
integer, intent(in) :: &
ip, &
el
real(pReal), intent(out) :: Tdot
end subroutine thermal_conduction_getSource
module function damage_nonlocal_getMobility(ce) result(M)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: M real(pReal), dimension(3,3) :: K
end function damage_nonlocal_getMobility end function homogenization_K_phi
module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce) module function homogenization_f_phi(phi,ce) result(f)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: & real(pReal), intent(in) :: phi
phi real(pReal) :: f
real(pReal), intent(out) :: & end function homogenization_f_phi
phiDot
end subroutine damage_nonlocal_getSourceAndItsTangent
module subroutine homogenization_set_phi(phi,ce) module subroutine homogenization_set_phi(phi,ce)
integer, intent(in) :: ce integer, intent(in) :: ce
@ -187,14 +178,14 @@ module homogenization
public :: & public :: &
homogenization_init, & homogenization_init, &
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
homogenization_thermal_mu_T, & homogenization_mu_T, &
thermal_conduction_getConductivity, & homogenization_K_T, &
thermal_conduction_getSource, & homogenization_f_T, &
damage_nonlocal_getMobility, &
damage_nonlocal_getSourceAndItsTangent, &
homogenization_set_phi, &
homogenization_thermal_setfield, & homogenization_thermal_setfield, &
homogenization_thermal_T, & homogenization_mu_phi, &
homogenization_K_phi, &
homogenization_f_phi, &
homogenization_set_phi, &
homogenization_forward, & homogenization_forward, &
homogenization_results, & homogenization_results, &
homogenization_restartRead, & homogenization_restartRead, &
@ -202,9 +193,6 @@ module homogenization
THERMAL_CONDUCTION_ID, & THERMAL_CONDUCTION_ID, &
DAMAGE_NONLOCAL_ID DAMAGE_NONLOCAL_ID
public :: &
damage_nonlocal_getDiffusion
contains contains
@ -314,7 +302,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
terminallyIll = .true. ! ...and kills all others terminallyIll = .true. ! ...and kills all others
endif endif
enddo enddo
call thermal_homogenize(ip,el)
enddo enddo
enddo enddo
!$OMP END DO !$OMP END DO
@ -447,32 +434,6 @@ subroutine homogenization_restartRead(fileHandle)
end subroutine homogenization_restartRead end subroutine homogenization_restartRead
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized non local damage diffusion tensor in reference configuration
!--------------------------------------------------------------------------------------------------
function damage_nonlocal_getDiffusion(ce)
integer, intent(in) :: ce
real(pReal), dimension(3,3) :: &
damage_nonlocal_getDiffusion
integer :: &
ho, &
co
ho = material_homogenizationID(ce)
damage_nonlocal_getDiffusion = 0.0_pReal
do co = 1, homogenization_Nconstituents(ho)
damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseID(co,ce)))
enddo
damage_nonlocal_getDiffusion = &
num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(ho),pReal)
end function damage_nonlocal_getDiffusion
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration !> @brief parses the homogenization part from the material configuration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

86
src/homogenization_damage.f90
View File

@ -26,8 +26,8 @@ submodule(homogenization) damage
type(tparameters), dimension(:), allocatable :: & type(tparameters), dimension(:), allocatable :: &
param param
contains
contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Allocate variables and set parameters. !> @brief Allocate variables and set parameters.
@ -38,14 +38,12 @@ module subroutine damage_init()
configHomogenizations, & configHomogenizations, &
configHomogenization, & configHomogenization, &
configHomogenizationDamage, & configHomogenizationDamage, &
num_generic, & num_generic
material_homogenization integer :: ho,Nmaterialpoints
integer :: ho
integer :: Ninstances,Nmaterialpoints,h
print'(/,a)', ' <<<+- homogenization:damage init -+>>>' print'(/,a)', ' <<<+- homogenization:damage init -+>>>'
print'(/,a)', ' <<<+- homogenization:damage:pass init -+>>>'
configHomogenizations => config_material%get('homogenization') configHomogenizations => config_material%get('homogenization')
allocate(param(configHomogenizations%length)) allocate(param(configHomogenizations%length))
@ -62,6 +60,10 @@ module subroutine damage_init()
#else #else
prm%output = configHomogenizationDamage%get_as1dString('output',defaultVal=emptyStringArray) prm%output = configHomogenizationDamage%get_as1dString('output',defaultVal=emptyStringArray)
#endif #endif
Nmaterialpoints = count(material_homogenizationAt == ho)
damageState_h(ho)%sizeState = 1
allocate(damageState_h(ho)%state0(1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState_h(ho)%state (1,Nmaterialpoints), source=1.0_pReal)
else else
prm%output = emptyStringArray prm%output = emptyStringArray
endif endif
@ -73,18 +75,7 @@ module subroutine damage_init()
num_generic => config_numerics%get('generic',defaultVal= emptyDict) num_generic => config_numerics%get('generic',defaultVal= emptyDict)
num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
Ninstances = count(damage_type == DAMAGE_nonlocal_ID) call pass_init()
material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
Nmaterialpoints = count(material_homogenizationAt == h)
damageState_h(h)%sizeState = 1
allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal)
enddo
end subroutine damage_init end subroutine damage_init
@ -100,53 +91,70 @@ module subroutine damage_partition(ce)
if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return
phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce)) phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce))
call phase_damage_set_phi(phi,1,ce) call phase_set_phi(phi,1,ce)
end subroutine damage_partition end subroutine damage_partition
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Returns homogenized nonlocal damage mobility !> @brief Homogenized damage viscosity.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function damage_nonlocal_getMobility(ce) result(M) module function homogenization_mu_phi(ce) result(mu)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: M real(pReal) :: mu
M = lattice_M(material_phaseID(1,ce))
end function damage_nonlocal_getMobility mu = phase_mu_phi(1,ce)
end function homogenization_mu_phi
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates homogenized damage driving forces !> @brief Homogenized damage conductivity/diffusivity in reference configuration.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce) module function homogenization_K_phi(ce) result(K)
integer, intent(in) :: ce
real(pReal), dimension(3,3) :: K
K = phase_K_phi(1,ce) &
* num_damage%charLength**2
end function homogenization_K_phi
!--------------------------------------------------------------------------------------------------
!> @brief Homogenized damage driving force.
!--------------------------------------------------------------------------------------------------
module function homogenization_f_phi(phi,ce) result(f)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
real(pReal), intent(out) :: & real(pReal) :: f
phiDot
phiDot = phase_damage_phi_dot(phi, 1, ce)
end subroutine damage_nonlocal_getSourceAndItsTangent f = phase_f_phi(phi, 1, ce)
end function homogenization_f_phi
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief updated nonlocal damage field with solution from damage phase field PDE !> @brief Set damage field.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine homogenization_set_phi(phi,ce) module subroutine homogenization_set_phi(phi,ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer :: & integer :: &
ho, & ho, &
en en
ho = material_homogenizationID(ce) ho = material_homogenizationID(ce)
en = material_homogenizationEntry(ce) en = material_homogenizationEntry(ce)
damagestate_h(ho)%state(1,en) = phi damagestate_h(ho)%state(1,en) = phi
@ -166,13 +174,13 @@ module subroutine damage_results(ho,group)
integer :: o integer :: o
associate(prm => param(ho)) associate(prm => param(ho))
outputsLoop: do o = 1,size(prm%output) outputsLoop: do o = 1,size(prm%output)
select case(prm%output(o)) select case(prm%output(o))
case ('phi') case ('phi')
call results_writeDataset(group,damagestate_h(ho)%state(1,:),prm%output(o),& call results_writeDataset(group,damagestate_h(ho)%state(1,:),prm%output(o),&
'damage indicator','-') 'damage indicator','-')
end select end select
enddo outputsLoop enddo outputsLoop
end associate end associate
end subroutine damage_results end subroutine damage_results

5
src/homogenization_damage_pass.f90
View File

@ -6,8 +6,9 @@ submodule(homogenization:damage) damage_pass
contains contains
module subroutine pass_init module subroutine pass_init()
print'(/,a)', ' <<<+- homogenization:damage:pass init -+>>>'
end subroutine pass_init end subroutine pass_init

61
src/homogenization_mechanical.f90
View File

@ -32,25 +32,6 @@ submodule(homogenization) mechanical
end subroutine RGC_partitionDeformation end subroutine RGC_partitionDeformation
module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: ho
end subroutine isostrain_averageStressAndItsTangent
module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: ho
end subroutine RGC_averageStressAndItsTangent
module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: & real(pReal), dimension(:,:,:), intent(in) :: &
@ -148,39 +129,21 @@ module subroutine mechanical_homogenize(dt,ce)
integer, intent(in) :: ce integer, intent(in) :: ce
integer :: co integer :: co
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce))) 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(dt,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)
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(dt,co,ce)
enddo
case (HOMOGENIZATION_NONE_ID) chosenHomogenization homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) &
homogenization_P(1:3,1:3,ce) = phase_P(1,ce) / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) 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)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
Ps(:,:,co) = phase_P(co,ce)
enddo
call isostrain_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, &
material_homogenizationID(ce))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
Ps(:,:,co) = phase_P(co,ce)
enddo
call RGC_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, &
material_homogenizationID(ce))
end select chosenHomogenization
end subroutine mechanical_homogenize end subroutine mechanical_homogenize

25
src/homogenization_mechanical_RGC.f90
View File

@ -89,7 +89,8 @@ module subroutine RGC_init(num_homogMech)
print'(/,a)', ' <<<+- homogenization:mechanical:RGC init -+>>>' print'(/,a)', ' <<<+- homogenization:mechanical:RGC init -+>>>'
print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_RGC_ID); flush(IO_STDOUT) print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_RGC_ID)
flush(IO_STDOUT)
print*, 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' print*, 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
@ -207,7 +208,7 @@ module subroutine RGC_partitionDeformation(F,avgF,ce)
integer :: iGrain,iFace,i,j,ho,en integer :: iGrain,iFace,i,j,ho,en
associate(prm => param(material_homogenizationID(ce))) associate(prm => param(material_homogenizationID(ce)))
ho = material_homogenizationID(ce) ho = material_homogenizationID(ce)
en = material_homogenizationEntry(ce) en = material_homogenizationEntry(ce)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -700,24 +701,6 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
end function RGC_updateState end function RGC_updateState
!--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities
!--------------------------------------------------------------------------------------------------
module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: ho
avgP = sum(P,3) /real(product(param(ho)%N_constituents),pReal)
dAvgPdAvgF = sum(dPdF,5)/real(product(param(ho)%N_constituents),pReal)
end subroutine RGC_averageStressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -802,7 +785,7 @@ pure function interfaceNormal(intFace,ho,en)
interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis
interfaceNormal = matmul(dst%orientation(1:3,1:3,en),interfaceNormal) ! map the normal vector into sample coordinate system (basis) interfaceNormal = matmul(dst%orientation(1:3,1:3,en),interfaceNormal) ! map the normal vector into sample coordinate system (basis)
end associate end associate
end function interfaceNormal end function interfaceNormal

83
src/homogenization_mechanical_isostrain.f90
View File

@ -6,21 +6,6 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization:mechanical) isostrain submodule(homogenization:mechanical) isostrain
enum, bind(c); enumerator :: &
parallel_ID, &
average_ID
end enum
type :: tParameters !< container type for internal constitutive parameters
integer :: &
N_constituents
integer(kind(average_ID)) :: &
mapping
end type
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -29,42 +14,21 @@ contains
module subroutine isostrain_init module subroutine isostrain_init
integer :: & integer :: &
h, & ho, &
Nmaterialpoints Nmaterialpoints
class(tNode), pointer :: &
material_homogenization, &
homog, &
homogMech
print'(/,a)', ' <<<+- homogenization:mechanical:isostrain init -+>>>' print'(/,a)', ' <<<+- homogenization:mechanical:isostrain init -+>>>'
print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID); flush(IO_STDOUT) print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
flush(IO_STDOUT)
material_homogenization => config_material%get('homogenization') do ho = 1, size(homogenization_type)
allocate(param(material_homogenization%length)) ! one container of parameters per homog if (homogenization_type(ho) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle
do h = 1, size(homogenization_type) Nmaterialpoints = count(material_homogenizationAt == ho)
if (homogenization_type(h) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle homogState(ho)%sizeState = 0
homog => material_homogenization%get(h) allocate(homogState(ho)%state0(0,Nmaterialpoints))
homogMech => homog%get('mechanical') allocate(homogState(ho)%state (0,Nmaterialpoints))
associate(prm => param(h))
prm%N_constituents = homogenization_Nconstituents(h)
select case(homogMech%get_asString('mapping',defaultVal = 'sum'))
case ('sum')
prm%mapping = parallel_ID
case ('avg')
prm%mapping = average_ID
case default
call IO_error(211,ext_msg='sum'//' (isostrain)')
end select
Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0
allocate(homogState(h)%state0 (0,Nmaterialpoints))
allocate(homogState(h)%state (0,Nmaterialpoints))
end associate
enddo enddo
@ -80,36 +44,9 @@ module subroutine isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
F = spread(avgF,3,size(F,3)) F = spread(avgF,3,size(F,3))
end subroutine isostrain_partitionDeformation end subroutine isostrain_partitionDeformation
!--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities
!--------------------------------------------------------------------------------------------------
module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: ho
associate(prm => param(ho))
select case (prm%mapping)
case (parallel_ID)
avgP = sum(P,3)
dAvgPdAvgF = sum(dPdF,5)
case (average_ID)
avgP = sum(P,3) /real(prm%N_constituents,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(prm%N_constituents,pReal)
end select
end associate
end subroutine isostrain_averageStressAndItsTangent
end submodule isostrain end submodule isostrain

21
src/homogenization_mechanical_pass.f90
View File

@ -14,25 +14,24 @@ contains
module subroutine pass_init module subroutine pass_init
integer :: & integer :: &
Ninstances, & ho, &
h, &
Nmaterialpoints Nmaterialpoints
print'(/,a)', ' <<<+- homogenization:mechanical:pass init -+>>>' print'(/,a)', ' <<<+- homogenization:mechanical:pass init -+>>>'
Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID) print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_NONE_ID)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) flush(IO_STDOUT)
do h = 1, size(homogenization_type) do ho = 1, size(homogenization_type)
if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle if(homogenization_type(ho) /= HOMOGENIZATION_NONE_ID) cycle
if(homogenization_Nconstituents(h) /= 1) & if(homogenization_Nconstituents(ho) /= 1) &
call IO_error(211,ext_msg='N_constituents (pass)') call IO_error(211,ext_msg='N_constituents (pass)')
Nmaterialpoints = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == ho)
homogState(h)%sizeState = 0 homogState(ho)%sizeState = 0
allocate(homogState(h)%state0 (0,Nmaterialpoints)) allocate(homogState(ho)%state0(0,Nmaterialpoints))
allocate(homogState(h)%state (0,Nmaterialpoints)) allocate(homogState(ho)%state (0,Nmaterialpoints))
enddo enddo

127
src/homogenization_thermal.f90
View File

@ -45,8 +45,6 @@ module subroutine thermal_init()
print'(/,a)', ' <<<+- homogenization:thermal init -+>>>' print'(/,a)', ' <<<+- homogenization:thermal init -+>>>'
print'(/,a)', ' <<<+- homogenization:thermal:pass init -+>>>'
configHomogenizations => config_material%get('homogenization') configHomogenizations => config_material%get('homogenization')
@ -71,6 +69,8 @@ module subroutine thermal_init()
end associate end associate
enddo enddo
call pass_init()
end subroutine thermal_init end subroutine thermal_init
@ -95,94 +95,70 @@ end subroutine thermal_partition
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Homogenize temperature rates !> @brief Homogenized thermal viscosity.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine thermal_homogenize(ip,el) module function homogenization_mu_T(ce) result(mu)
integer, intent(in) :: ip,el integer, intent(in) :: ce
real(pReal) :: mu
!call phase_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el) integer :: co
end subroutine thermal_homogenize
mu = phase_mu_T(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
mu = mu + phase_mu_T(co,ce)
enddo
mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function homogenization_mu_T
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return homogenized thermal conductivity in reference configuration !> @brief Homogenized thermal conductivity in reference configuration.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function thermal_conduction_getConductivity(ce) result(K) module function homogenization_K_T(ce) result(K)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), dimension(3,3) :: K real(pReal), dimension(3,3) :: K
integer :: & integer :: co
co
K = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) K = phase_K_T(1,ce)
K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseID(co,ce))) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
K = K + phase_K_T(co,ce)
enddo enddo
K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function thermal_conduction_getConductivity end function homogenization_K_T
module function homogenization_thermal_mu_T(ce) result(mu_T)
integer, intent(in) :: ce
real(pReal) :: mu_T
mu_T = c_P(ce) * rho(ce)
end function homogenization_thermal_mu_T
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized specific heat capacity !> @brief Homogenized heat generation rate.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function c_P(ce) module function homogenization_f_T(ce) result(f)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: c_P real(pReal) :: f
integer :: co integer :: co
c_P = lattice_c_p(material_phaseID(1,ce)) f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
c_P = c_P + lattice_c_p(material_phaseID(co,ce)) f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
enddo enddo
c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function c_P end function homogenization_f_T
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized mass density !> @brief Set thermal field and its rate (T and dot_T).
!--------------------------------------------------------------------------------------------------
function rho(ce)
integer, intent(in) :: ce
real(pReal) :: rho
integer :: co
rho = lattice_rho(material_phaseID(1,ce))
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
rho = rho + lattice_rho(material_phaseID(co,ce))
enddo
rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function rho
!--------------------------------------------------------------------------------------------------
!> @brief Set thermal field and its rate (T and dot_T)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine homogenization_thermal_setField(T,dot_T, ce) module subroutine homogenization_thermal_setField(T,dot_T, ce)
@ -197,7 +173,6 @@ module subroutine homogenization_thermal_setField(T,dot_T, ce)
end subroutine homogenization_thermal_setField end subroutine homogenization_thermal_setField
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -219,44 +194,4 @@ module subroutine thermal_results(ho,group)
end subroutine thermal_results end subroutine thermal_results
module function homogenization_thermal_T(ce) result(T)
integer, intent(in) :: ce
real(pReal) :: T
T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce))
end function homogenization_thermal_T
!--------------------------------------------------------------------------------------------------
!> @brief return heat generation rate
!--------------------------------------------------------------------------------------------------
module subroutine thermal_conduction_getSource(Tdot, ip, el)
integer, intent(in) :: &
ip, &
el
real(pReal), intent(out) :: &
Tdot
integer :: co, ho,ph,me
real(pReal) :: dot_T_temp
ho = material_homogenizationAt(el)
Tdot = 0.0_pReal
do co = 1, homogenization_Nconstituents(ho)
ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el)
call phase_thermal_getRate(dot_T_temp, ph,me)
Tdot = Tdot + dot_T_temp
enddo
Tdot = Tdot/real(homogenization_Nconstituents(ho),pReal)
end subroutine thermal_conduction_getSource
end submodule thermal end submodule thermal

5
src/homogenization_thermal_pass.f90
View File

@ -6,8 +6,9 @@ submodule(homogenization:thermal) thermal_pass
contains contains
module subroutine pass_init module subroutine pass_init()
print'(/,a)', ' <<<+- homogenization:thermal:pass init -+>>>'
end subroutine pass_init end subroutine pass_init

45
src/lattice.f90
View File

@ -394,13 +394,11 @@ module lattice
! SHOULD NOT BE PART OF LATTICE BEGIN ! SHOULD NOT BE PART OF LATTICE BEGIN
real(pReal), dimension(:), allocatable, public, protected :: & real(pReal), dimension(:), allocatable, public, protected :: &
lattice_mu, lattice_nu, & lattice_mu, lattice_nu, &
lattice_M, & lattice_mu_phi, &
lattice_rho, & lattice_rho, &
lattice_c_p lattice_c_p
real(pReal), dimension(:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_C66, & lattice_C66
lattice_K, &
lattice_D
integer(kind(lattice_UNDEFINED_ID)), dimension(:), allocatable, public, protected :: & integer(kind(lattice_UNDEFINED_ID)), dimension(:), allocatable, public, protected :: &
lattice_structure lattice_structure
! SHOULD NOT BE PART OF LATTICE END ! SHOULD NOT BE PART OF LATTICE END
@ -458,23 +456,19 @@ subroutine lattice_init
phases, & phases, &
phase, & phase, &
mech, & mech, &
elasticity, & elasticity
thermal, &
damage
print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT)
! SHOULD NOT BE PART OF LATTICE BEGIN
phases => config_material%get('phase') phases => config_material%get('phase')
Nphases = phases%length Nphases = phases%length
allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID) allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID)
allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) allocate(lattice_C66(6,6,Nphases), source=0.0_pReal)
allocate(lattice_K (3,3,Nphases), source=0.0_pReal) allocate(lattice_rho, &
allocate(lattice_D (3,3,Nphases), source=0.0_pReal)
allocate(lattice_M,&
lattice_rho,lattice_c_p, &
lattice_mu, lattice_nu,& lattice_mu, lattice_nu,&
source=[(0.0_pReal,i=1,Nphases)]) source=[(0.0_pReal,i=1,Nphases)])
@ -522,32 +516,7 @@ subroutine lattice_init
enddo enddo
lattice_rho(ph) = phase%get_asFloat('rho', defaultVal=0.0_pReal) lattice_rho(ph) = phase%get_asFloat('rho', defaultVal=0.0_pReal)
! SHOULD NOT BE PART OF LATTICE END
! SHOULD NOT BE PART OF LATTICE BEGIN
if (phase%contains('thermal')) then
thermal => phase%get('thermal')
lattice_K(1,1,ph) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)
lattice_K(2,2,ph) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal)
lattice_K(3,3,ph) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)
lattice_K(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,ph), &
phase%get_asString('lattice'))
lattice_c_p(ph) = thermal%get_asFloat('c_p', defaultVal=0.0_pReal)
endif
if (phase%contains('damage')) then
damage => phase%get('damage')
damage => damage%get(1)
lattice_D(1,1,ph) = damage%get_asFloat('D_11',defaultVal=0.0_pReal)
lattice_D(2,2,ph) = damage%get_asFloat('D_22',defaultVal=0.0_pReal)
lattice_D(3,3,ph) = damage%get_asFloat('D_33',defaultVal=0.0_pReal)
lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), &
phase%get_asString('lattice'))
lattice_M(ph) = damage%get_asFloat('M',defaultVal=0.0_pReal)
endif
! SHOULD NOT BE PART OF LATTICE END
call selfTest call selfTest

68
src/phase.f90
View File

@ -145,26 +145,22 @@ module phase
real(pReal), dimension(3,3) :: L_p real(pReal), dimension(3,3) :: L_p
end function mechanical_L_p end function mechanical_L_p
module function phase_F(co,ce) result(F)
integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: F
end function phase_F
module function mechanical_F_e(ph,me) result(F_e) module function mechanical_F_e(ph,me) result(F_e)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: F_e real(pReal), dimension(3,3) :: F_e
end function mechanical_F_e end function mechanical_F_e
module function phase_F(co,ce) result(F)
integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: F
end function phase_F
module function phase_P(co,ce) result(P) module function phase_P(co,ce) result(P)
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: P real(pReal), dimension(3,3) :: P
end function phase_P end function phase_P
module function phase_damage_get_phi(co,ip,el) result(phi)
integer, intent(in) :: co, ip, el
real(pReal) :: phi
end function phase_damage_get_phi
module function thermal_T(ph,me) result(T) module function thermal_T(ph,me) result(T)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal) :: T real(pReal) :: T
@ -188,13 +184,35 @@ module phase
module subroutine phase_thermal_setField(T,dot_T, co,ce) module subroutine phase_thermal_setField(T,dot_T, co,ce)
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
integer, intent(in) :: ce, co integer, intent(in) :: co, ce
end subroutine phase_thermal_setField end subroutine phase_thermal_setField
module subroutine phase_damage_set_phi(phi,co,ce) module subroutine phase_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi real(pReal), intent(in) :: phi
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
end subroutine phase_damage_set_phi end subroutine phase_set_phi
module function phase_mu_phi(co,ce) result(mu)
integer, intent(in) :: co, ce
real(pReal) :: mu
end function phase_mu_phi
module function phase_K_phi(co,ce) result(K)
integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: K
end function phase_K_phi
module function phase_mu_T(co,ce) result(mu)
integer, intent(in) :: co, ce
real(pReal) :: mu
end function phase_mu_T
module function phase_K_T(co,ce) result(K)
integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: K
end function phase_K_T
! == cleaned:end =================================================================================== ! == cleaned:end ===================================================================================
@ -227,19 +245,18 @@ module phase
end function phase_homogenizedC end function phase_homogenizedC
module function phase_damage_phi_dot(phi,co,ce) result(phi_dot) module function phase_f_phi(phi,co,ce) result(f)
integer, intent(in) :: ce,co integer, intent(in) :: ce,co
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi !< damage parameter phi !< damage parameter
real(pReal) :: & real(pReal) :: &
phi_dot f
end function phase_damage_phi_dot end function phase_f_phi
module subroutine phase_thermal_getRate(TDot, ph,me) module function phase_f_T(ph,me) result(f)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), intent(out) :: & real(pReal) :: f
TDot end function phase_f_T
end subroutine phase_thermal_getRate
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
integer, intent(in) :: & integer, intent(in) :: &
@ -300,8 +317,12 @@ module phase
public :: & public :: &
phase_init, & phase_init, &
phase_homogenizedC, & phase_homogenizedC, &
phase_damage_phi_dot, & phase_f_phi, &
phase_thermal_getRate, & phase_f_T, &
phase_K_phi, &
phase_K_T, &
phase_mu_phi, &
phase_mu_T, &
phase_results, & phase_results, &
phase_allocateState, & phase_allocateState, &
phase_forward, & phase_forward, &
@ -318,8 +339,7 @@ module phase
phase_restartRead, & phase_restartRead, &
integrateDamageState, & integrateDamageState, &
phase_thermal_setField, & phase_thermal_setField, &
phase_damage_set_phi, & phase_set_phi, &
phase_damage_get_phi, &
phase_P, & phase_P, &
phase_set_F, & phase_set_F, &
phase_F phase_F

135
src/phase_damage.f90
View File

@ -2,6 +2,12 @@
!> @brief internal microstructure state for all damage sources and kinematics constitutive models !> @brief internal microstructure state for all damage sources and kinematics constitutive models
!---------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------
submodule(phase) damage submodule(phase) damage
type :: tDamageParameters
real(pReal) :: mu = 0.0_pReal !< viscosity
real(pReal), dimension(3,3) :: K = 0.0_pReal !< conductivity/diffusivity
end type tDamageParameters
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
DAMAGE_UNDEFINED_ID, & DAMAGE_UNDEFINED_ID, &
DAMAGE_ISOBRITTLE_ID, & DAMAGE_ISOBRITTLE_ID, &
@ -16,13 +22,12 @@ submodule(phase) damage
end type tDataContainer end type tDataContainer
integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: &
phase_source !< active sources mechanisms of each phase phase_source !< active sources mechanisms of each phase
integer, dimension(:), allocatable :: &
phase_Nsources
type(tDataContainer), dimension(:), allocatable :: current type(tDataContainer), dimension(:), allocatable :: current
type(tDamageParameters), dimension(:), allocatable :: param
interface interface
module function anisobrittle_init() result(mySources) module function anisobrittle_init() result(mySources)
@ -100,18 +105,19 @@ module subroutine damage_init
class(tNode), pointer :: & class(tNode), pointer :: &
phases, & phases, &
phase, & phase, &
sources sources, &
source
logical:: damage_active
print'(/,a)', ' <<<+- phase:damage init -+>>>' print'(/,a)', ' <<<+- phase:damage init -+>>>'
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(current(phases%length)) allocate(current(phases%length))
allocate(damageState (phases%length)) allocate(damageState (phases%length))
allocate(phase_Nsources(phases%length),source = 0) allocate(param(phases%length))
damage_active = .false.
do ph = 1,phases%length do ph = 1,phases%length
Nmembers = count(material_phaseID == ph) Nmembers = count(material_phaseID == ph)
@ -122,14 +128,21 @@ module subroutine damage_init
phase => phases%get(ph) phase => phases%get(ph)
sources => phase%get('damage',defaultVal=emptyList) sources => phase%get('damage',defaultVal=emptyList)
if (sources%length > 1) error stop if (sources%length > 1) error stop
phase_Nsources(ph) = sources%length if (sources%length == 1) then
damage_active = .true.
source => sources%get(1)
param(ph)%mu = source%get_asFloat('M',defaultVal=0.0_pReal)
param(ph)%K(1,1) = source%get_asFloat('D_11',defaultVal=0.0_pReal)
param(ph)%K(2,2) = source%get_asFloat('D_22',defaultVal=0.0_pReal)
param(ph)%K(3,3) = source%get_asFloat('D_33',defaultVal=0.0_pReal)
param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase%get_asString('lattice'))
endif
enddo enddo
allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID) allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID)
! initialize source mechanisms if (damage_active) then
if(maxval(phase_Nsources) /= 0) then
where(isobrittle_init() ) phase_source = DAMAGE_ISOBRITTLE_ID where(isobrittle_init() ) phase_source = DAMAGE_ISOBRITTLE_ID
where(isoductile_init() ) phase_source = DAMAGE_ISODUCTILE_ID where(isoductile_init() ) phase_source = DAMAGE_ISODUCTILE_ID
where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID
@ -142,13 +155,13 @@ end subroutine damage_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief returns local part of nonlocal damage driving force !< @brief returns local part of nonlocal damage driving force
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function phase_damage_phi_dot(phi,co,ce) result(phi_dot) module function phase_f_phi(phi,co,ce) result(f)
integer, intent(in) :: ce,co integer, intent(in) :: ce,co
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi !< damage parameter phi !< damage parameter
real(pReal) :: & real(pReal) :: &
phi_dot f
integer :: & integer :: &
ph, & ph, &
@ -159,13 +172,13 @@ module function phase_damage_phi_dot(phi,co,ce) result(phi_dot)
select case(phase_source(ph)) select case(phase_source(ph))
case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID,DAMAGE_ANISODUCTILE_ID) case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID,DAMAGE_ANISODUCTILE_ID)
phi_dot = 1.0_pReal & f = 1.0_pReal &
- phi*damageState(ph)%state(1,en) - phi*damageState(ph)%state(1,en)
case default case default
phi_dot = 0.0_pReal f = 0.0_pReal
end select end select
end function phase_damage_phi_dot end function phase_f_phi
@ -276,30 +289,25 @@ module subroutine damage_results(group,ph)
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer, intent(in) :: ph integer, intent(in) :: ph
integer :: so
sourceLoop: do so = 1, phase_Nsources(ph)
if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) & if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) &
call results_closeGroup(results_addGroup(group//'damage')) call results_closeGroup(results_addGroup(group//'damage'))
sourceType: select case (phase_source(ph)) sourceType: select case (phase_source(ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE_ID) sourceType
call isobrittle_results(ph,group//'damage/') call isobrittle_results(ph,group//'damage/')
case (DAMAGE_ISODUCTILE_ID) sourceType case (DAMAGE_ISODUCTILE_ID) sourceType
call isoductile_results(ph,group//'damage/') call isoductile_results(ph,group//'damage/')
case (DAMAGE_ANISOBRITTLE_ID) sourceType case (DAMAGE_ANISOBRITTLE_ID) sourceType
call anisobrittle_results(ph,group//'damage/') call anisobrittle_results(ph,group//'damage/')
case (DAMAGE_ANISODUCTILE_ID) sourceType case (DAMAGE_ANISODUCTILE_ID) sourceType
call anisoductile_results(ph,group//'damage/') call anisoductile_results(ph,group//'damage/')
end select sourceType end select sourceType
enddo SourceLoop
end subroutine damage_results end subroutine damage_results
@ -339,6 +347,33 @@ function phase_damage_collectDotState(ph,me) result(broken)
end function phase_damage_collectDotState end function phase_damage_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief Damage viscosity.
!--------------------------------------------------------------------------------------------------
module function phase_mu_phi(co,ce) result(mu)
integer, intent(in) :: co, ce
real(pReal) :: mu
mu = param(material_phaseID(co,ce))%mu
end function phase_mu_phi
!--------------------------------------------------------------------------------------------------
!> @brief Damage conductivity/diffusivity in reference configuration.
!--------------------------------------------------------------------------------------------------
module function phase_K_phi(co,ce) result(K)
integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: K
K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%K)
end function phase_K_phi
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
@ -362,19 +397,19 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
if (damageState(ph)%sizeState == 0) return if (damageState(ph)%sizeState == 0) return
sourceType: select case (phase_source(ph)) sourceType: select case (phase_source(ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE_ID) sourceType
call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me)))
if(.not. broken) then if(.not. broken) then
myOffset = damageState(ph)%offsetDeltaState myOffset = damageState(ph)%offsetDeltaState
mySize = damageState(ph)%sizeDeltaState mySize = damageState(ph)%sizeDeltaState
damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = & damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = &
damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me) damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me)
endif endif
end select sourceType end select sourceType
end function phase_damage_deltaState end function phase_damage_deltaState
@ -411,7 +446,7 @@ end function source_active
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Set damage parameter !< @brief Set damage parameter
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine phase_damage_set_phi(phi,co,ce) module subroutine phase_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi real(pReal), intent(in) :: phi
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
@ -419,17 +454,7 @@ module subroutine phase_damage_set_phi(phi,co,ce)
current(material_phaseID(co,ce))%phi(material_phaseEntry(co,ce)) = phi current(material_phaseID(co,ce))%phi(material_phaseEntry(co,ce)) = phi
end subroutine phase_damage_set_phi end subroutine phase_set_phi
module function phase_damage_get_phi(co,ip,el) result(phi)
integer, intent(in) :: co, ip, el
real(pReal) :: phi
phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el))
end function phase_damage_get_phi
module function damage_phi(ph,me) result(phi) module function damage_phi(ph,me) result(phi)

87
src/phase_thermal.f90
View File

@ -3,6 +3,11 @@
!---------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------
submodule(phase) thermal submodule(phase) thermal
type :: tThermalParameters
real(pReal) :: C_p = 0.0_pReal !< heat capacity
real(pReal), dimension(3,3) :: K = 0.0_pReal !< thermal conductivity
end type tThermalParameters
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
thermal_Nsources thermal_Nsources
@ -21,7 +26,9 @@ submodule(phase) thermal
integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: &
thermal_source thermal_source
type(tDataContainer), dimension(:), allocatable :: current ! ?? not very telling name. Better: "field" ?? type(tDataContainer), dimension(:), allocatable :: current ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(me) reads quite good
type(tThermalParameters), dimension(:), allocatable :: param
integer :: thermal_source_maxSizeDotState integer :: thermal_source_maxSizeDotState
@ -45,21 +52,19 @@ submodule(phase) thermal
me me
end subroutine externalheat_dotState end subroutine externalheat_dotState
module subroutine dissipation_getRate(TDot, ph,me) module function dissipation_f_T(ph,me) result(f_T)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me me
real(pReal), intent(out) :: & real(pReal) :: f_T
TDot end function dissipation_f_T
end subroutine dissipation_getRate
module subroutine externalheat_getRate(TDot, ph,me) module function externalheat_f_T(ph,me) result(f_T)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me me
real(pReal), intent(out) :: & real(pReal) :: f_T
TDot end function externalheat_f_T
end subroutine externalheat_getRate
end interface end interface
@ -87,6 +92,7 @@ module subroutine thermal_init(phases)
allocate(thermalState(phases%length)) allocate(thermalState(phases%length))
allocate(thermal_Nsources(phases%length),source = 0) allocate(thermal_Nsources(phases%length),source = 0)
allocate(param(phases%length))
do ph = 1, phases%length do ph = 1, phases%length
Nmembers = count(material_phaseID == ph) Nmembers = count(material_phaseID == ph)
@ -94,9 +100,17 @@ module subroutine thermal_init(phases)
allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal) allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
phase => phases%get(ph) phase => phases%get(ph)
thermal => phase%get('thermal',defaultVal=emptyDict) thermal => phase%get('thermal',defaultVal=emptyDict)
param(ph)%C_p = thermal%get_asFloat('c_p',defaultVal=0.0_pReal)
if (param(ph)%C_p <= 0) param(ph)%C_p = thermal%get_asFloat('C_p',defaultVal=0.0_pReal)
param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)
param(ph)%K(2,2) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal)
param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)
param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase%get_asString('lattice'))
sources => thermal%get('source',defaultVal=emptyList) sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length thermal_Nsources(ph) = sources%length
allocate(thermalstate(ph)%p(thermal_Nsources(ph))) allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
enddo enddo
allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
@ -123,35 +137,31 @@ end subroutine thermal_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief calculates thermal dissipation rate !< @brief calculates thermal dissipation rate
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine phase_thermal_getRate(TDot, ph,me) module function phase_f_T(ph,me) result(f)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), intent(out) :: & real(pReal) :: f
TDot
real(pReal) :: &
my_Tdot
integer :: &
so
TDot = 0.0_pReal integer :: so
f = 0.0_pReal
do so = 1, thermal_Nsources(ph) do so = 1, thermal_Nsources(ph)
select case(thermal_source(so,ph)) select case(thermal_source(so,ph))
case (THERMAL_DISSIPATION_ID) case (THERMAL_DISSIPATION_ID)
call dissipation_getRate(my_Tdot, ph,me) f = f + dissipation_f_T(ph,me)
case (THERMAL_EXTERNALHEAT_ID) case (THERMAL_EXTERNALHEAT_ID)
call externalheat_getRate(my_Tdot, ph,me) f = f + externalheat_f_T(ph,me)
case default
my_Tdot = 0.0_pReal
end select end select
Tdot = Tdot + my_Tdot
enddo enddo
end subroutine phase_thermal_getRate end function phase_f_T
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -179,6 +189,35 @@ function phase_thermal_collectDotState(ph,me) result(broken)
end function phase_thermal_collectDotState end function phase_thermal_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief Damage viscosity.
!--------------------------------------------------------------------------------------------------
module function phase_mu_T(co,ce) result(mu)
integer, intent(in) :: co, ce
real(pReal) :: mu
mu = lattice_rho(material_phaseID(co,ce)) &
* param(material_phaseID(co,ce))%C_p
end function phase_mu_T
!--------------------------------------------------------------------------------------------------
!> @brief Damage conductivity/diffusivity in reference configuration.
!--------------------------------------------------------------------------------------------------
module function phase_K_T(co,ce) result(K)
integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: K
K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%K)
end function phase_K_T
module function thermal_stress(Delta_t,ph,me) result(converged_) ! ?? why is this called "stress" when it seems closer to "updateState" ?? module function thermal_stress(Delta_t,ph,me) result(converged_) ! ?? why is this called "stress" when it seems closer to "updateState" ??
real(pReal), intent(in) :: Delta_t real(pReal), intent(in) :: Delta_t

10
src/phase_thermal_dissipation.f90
View File

@ -69,17 +69,17 @@ end function dissipation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Ninstancess dissipation rate !> @brief Ninstancess dissipation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine dissipation_getRate(TDot, ph,me) module function dissipation_f_T(ph,me) result(f_T)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), intent(out) :: & real(pReal) :: &
TDot f_T
associate(prm => param(ph)) associate(prm => param(ph))
TDot = prm%kappa*sum(abs(mechanical_S(ph,me)*mechanical_L_p(ph,me))) f_T = prm%kappa*sum(abs(mechanical_S(ph,me)*mechanical_L_p(ph,me)))
end associate end associate
end subroutine dissipation_getRate end function dissipation_f_T
end submodule dissipation end submodule dissipation

12
src/phase_thermal_externalheat.f90
View File

@ -100,13 +100,13 @@ end subroutine externalheat_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate !> @brief returns local heat generation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine externalheat_getRate(TDot, ph, me) module function externalheat_f_T(ph,me) result(f_T)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me me
real(pReal), intent(out) :: & real(pReal) :: &
TDot f_T
integer :: & integer :: &
so, interval so, interval
@ -122,12 +122,12 @@ module subroutine externalheat_getRate(TDot, ph, me)
if ( (frac_time < 0.0_pReal .and. interval == 1) & if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + & f_T = prm%f_T(interval ) * (1.0_pReal - frac_time) + &
prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside me bounds ! ...or extrapolate if outside me bounds
enddo enddo
end associate end associate
end subroutine externalheat_getRate end function externalheat_f_T
end submodule externalheat end submodule externalheat