diff --git a/PRIVATE b/PRIVATE index 129812414..1490f9741 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 1298124143e7e2901d0b9c2e79ab6388cb78a1e3 +Subproject commit 1490f97417664d6ae11d7ceafb6b98c9cc2dade1 diff --git a/examples/.gitignore b/examples/.gitignore index 93d78295b..b832fb271 100644 --- a/examples/.gitignore +++ b/examples/.gitignore @@ -3,3 +3,6 @@ *.xdmf *.sta *.vt* +*.out +*.sts +*.t16 diff --git a/examples/Marc/material.config b/examples/Marc/material.config deleted file mode 100644 index 46ea44367..000000000 --- a/examples/Marc/material.config +++ /dev/null @@ -1,429 +0,0 @@ -#-------------------# - -#-------------------# - -{../ConfigFiles/Homogenization_None_Dummy.config} - -#-------------------# - -#-------------------# - -[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 - -#-------------------# - -#-------------------# - -[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 - - -#-------------------# - -#-------------------# - -{../ConfigFiles/Phase_Phenopowerlaw_Aluminum.config} - -{../ConfigFiles/Phase_Isotropic_AluminumIsotropic.config} diff --git a/examples/Marc/material.yaml b/examples/Marc/material.yaml new file mode 100644 index 000000000..5c893ac52 --- /dev/null +++ b/examples/Marc/material.yaml @@ -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] diff --git a/examples/Marc/rotation.mud b/examples/Marc/rotation.mud deleted file mode 100644 index 21eff7974..000000000 Binary files a/examples/Marc/rotation.mud and /dev/null differ diff --git a/examples/Marc/sheet_r-value.dat b/examples/Marc/sheet_r-value.dat new file mode 100644 index 000000000..529f4d31a --- /dev/null +++ b/examples/Marc/sheet_r-value.dat @@ -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 +$................... diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index 14e21144a..9fb07fc1c 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -1,7 +1,6 @@ import subprocess import shlex import re -import io import os from pathlib import Path @@ -41,13 +40,9 @@ class Marc: return path_tools - def submit_job(self, - model, - job = 'job1', - logfile = False, + def submit_job(self, model, job, compile = False, - optimization = '', - ): + optimization = ''): usersub = Path(os.environ['DAMASK_ROOT'])/'src/DAMASK_Marc' usersub = usersub.parent/(usersub.name + ('.f90' if compile else '.marc')) @@ -64,22 +59,16 @@ class Marc: ' -prog ' + str(usersub.with_suffix('')) print(cmd) - if logfile is not None: - 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) + ret = subprocess.run(shlex.split(cmd),capture_output=True) 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): + print(ret.stderr.decode()) + print(ret.stdout.decode()) raise RuntimeError('Marc simulation failed (unknown return value)') - if v != 3004: - raise RuntimeError(f'Marc simulation failed ({v})') diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 4961b99d4..d8f4a0c3c 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -150,15 +150,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS H 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 ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll - elCP = mesh_FEM2DAMASK_elem(elFE) - - ma = (elCP-1) * discretization_nIPs + ip + elCP = discretization_Marc_FEM2DAMASK_elem(elFE) + ce = discretization_Marc_FEM2DAMASK_cell(ip,elFE) if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then 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_inp !end select chosenThermal1 - homogenization_F0(1:3,1:3,ma) = ffn - homogenization_F(1:3,1:3,ma) = ffn1 + homogenization_F0(1:3,1:3,ce) = ffn + homogenization_F(1:3,1:3,ce) = ffn1 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 ! translate from P to sigma - Kirchhoff = matmul(homogenization_P(1:3,1:3,ma), transpose(homogenization_F(1:3,1:3,ma))) - J_inverse = 1.0_pReal / math_det33(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,ce)) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) ! translate from dP/dF to dCS/dE 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 H(i,j,k,l) = H(i,j,k,l) & - + homogenization_F(j,m,ma) * homogenization_F(l,n,ma) & - * homogenization_dPdF(i,m,k,n,ma) & - - math_delta(j,l) * homogenization_F(i,m,ma) * homogenization_P(k,m,ma) & + + homogenization_F(j,m,ce) * homogenization_F(l,n,ce) & + * homogenization_dPdF(i,m,k,n,ce) & + - 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) & + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k)) enddo; enddo; enddo; enddo; enddo; enddo diff --git a/src/DAMASK_Marc.f90 b/src/DAMASK_Marc.f90 index e417be2fa..910ca86c0 100644 --- a/src/DAMASK_Marc.f90 +++ b/src/DAMASK_Marc.f90 @@ -218,7 +218,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & logical :: cutBack real(pReal), dimension(6) :: stress 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(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 elseif (lovl == 6) then ! stress requested by marc computationMode = CPFEM_CALCRESULTS - cp_en = mesh_FEM2DAMASK_elem(m(1)) if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" terminallyIll = .false. 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) :: & f + f(1) = homogenization_f_T(discretization_Marc_FEM2DAMASK_cell(n(3),n(1))) f(2) = 0.0_pReal - call thermal_conduction_getSource(f(1), n(3),mesh_FEM2DAMASK_elem(n(1))) end subroutine flux @@ -370,11 +369,11 @@ subroutine uedinc(inc,incsub) if (inc > inc_written) then - allocate(d_n(3,count(mesh_FEM2DAMASK_node /= -1))) - do n = lbound(mesh_FEM2DAMASK_node,1), ubound(mesh_FEM2DAMASK_node,1) - if (mesh_FEM2DAMASK_node(n) /= -1) then - call nodvar(1,n,d_n(1:3,mesh_FEM2DAMASK_node(n)),nqncomp,nqdatatype) - if(nqncomp == 2) d_n(3,mesh_FEM2DAMASK_node(n)) = 0.0_pReal + allocate(d_n(3,count(discretization_Marc_FEM2DAMASK_node /= -1))) + do n = lbound(discretization_Marc_FEM2DAMASK_node,1), ubound(discretization_Marc_FEM2DAMASK_node,1) + if (discretization_Marc_FEM2DAMASK_node(n) /= -1) then + call nodvar(1,n,d_n(1:3,discretization_Marc_FEM2DAMASK_node(n)),nqncomp,nqdatatype) + if(nqncomp == 2) d_n(3,discretization_Marc_FEM2DAMASK_node(n)) = 0.0_pReal endif enddo diff --git a/src/Marc/discretization_Marc.f90 b/src/Marc/discretization_Marc.f90 index d92623215..c9f099758 100644 --- a/src/Marc/discretization_Marc.f90 +++ b/src/Marc/discretization_Marc.f90 @@ -24,8 +24,8 @@ module discretization_marc mesh_unitlength !< physical length of one unit in mesh MD: needs systematic_name integer, dimension(:), allocatable, public, protected :: & - mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID MD: Needs systematic name - mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID MD: needs systematic_name + discretization_Marc_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID + discretization_Marc_FEM2DAMASK_node !< DAMASK node ID for Marc node ID type tCellNodeDefinition @@ -39,8 +39,9 @@ module discretization_marc connectivity_cell !< cell connectivity for each element,ip/cell public :: & - discretization_marc_init, & - discretization_marc_UpdateNodeAndIpCoords + discretization_Marc_init, & + discretization_Marc_updateNodeAndIpCoords, & + discretization_Marc_FEM2DAMASK_cell contains @@ -48,7 +49,7 @@ contains !> @brief initializes the mesh by calling all necessary private routines the mesh module !! Order and routines strongly depend on type of solver !-------------------------------------------------------------------------------------------------- -subroutine discretization_marc_init +subroutine discretization_Marc_init real(pReal), dimension(:,:), allocatable :: & node0_elem, & !< node x,y,z coordinates (initially!) @@ -96,7 +97,7 @@ subroutine discretization_marc_init call buildCells(connectivity_cell,cellNodeDefinition,& elem,connectivity_elem) node0_cell = buildCellNodes(node0_elem) - + IP_reshaped = buildIPcoordinates(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_results -end subroutine discretization_marc_init +end subroutine discretization_Marc_init !-------------------------------------------------------------------------------------------------- !> @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(:,:), 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_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, & nElems,inputFile) - call inputRead_mapElems(mesh_FEM2DAMASK_elem,& + call inputRead_mapElems(discretization_Marc_FEM2DAMASK_elem,& nElems,elem%nNodes,inputFile) - call inputRead_mapNodes(mesh_FEM2DAMASK_node,& + call inputRead_mapNodes(discretization_Marc_FEM2DAMASK_node,& nNodes,inputFile) call inputRead_elemNodes(node0_elem, & @@ -532,7 +550,7 @@ subroutine inputRead_elemNodes(nodes, & if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then chunkPos = [4,1,10,11,30,31,50,51,70] 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 nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1) enddo @@ -653,11 +671,11 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent) j = 0 do i = 1,nElem 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 do k = 1,chunkPos(1)-2 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 nNodesAlreadyRead = chunkPos(1) - 2 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)) do k = 1,chunkPos(1) 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 nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) enddo @@ -718,7 +736,7 @@ subroutine inputRead_material(materialAt,& if (initialcondTableStyle == 2) m = m + 2 contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements do i = 1,contInts(1) - e = mesh_FEM2DAMASK_elem(contInts(1+i)) + e = discretization_Marc_FEM2DAMASK_elem(contInts(1+i)) materialAt(e) = myVal enddo if (initialcondTableStyle == 0) m = m + 1 diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index bb6fa9d8d..ad614d8c5 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -74,8 +74,10 @@ subroutine discretization_grid_init(restart) allocate(materialAt_global(0)) ! needed for IntelMPI endif + call MPI_Bcast(grid,3,MPI_INTEGER,0,PETSC_COMM_WORLD, ierr) 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) if (ierr /= 0) error stop 'MPI error' call MPI_Bcast(origin,3,MPI_DOUBLE,0,PETSC_COMM_WORLD, ierr) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 002a15708..967ddb73a 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -15,8 +15,8 @@ module grid_damage_spectral use IO use spectral_utilities use discretization_grid - use YAML_types use homogenization + use YAML_types use config implicit none @@ -61,7 +61,7 @@ contains !> @brief allocates all neccessary fields and fills them with data ! ToDo: Restart not implemented !-------------------------------------------------------------------------------------------------- -subroutine grid_damage_spectral_init +subroutine grid_damage_spectral_init() PetscInt, dimension(0:worldsize-1) :: localK DM :: damage_grid @@ -88,10 +88,10 @@ subroutine grid_damage_spectral_init num_generic => config_numerics%get('generic',defaultVal=emptyDict) 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%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_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_rtol') + 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%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') !-------------------------------------------------------------------------------------------------- ! 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_lastInc(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 updateReference @@ -216,21 +217,21 @@ end function grid_damage_spectral_solution subroutine grid_damage_spectral_forward(cutBack) logical, intent(in) :: cutBack - integer :: i, j, k, ce + integer :: i, j, k, ce DM :: dm_local - PetscScalar, dimension(:,:,:), pointer :: x_scal - PetscErrorCode :: ierr + PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscErrorCode :: ierr if (cutBack) then phi_current = phi_lastInc phi_stagInc = phi_lastInc !-------------------------------------------------------------------------------------------------- ! reverting damage field state - ce = 0 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 x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current 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) ce = ce + 1 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 PetscErrorCode :: ierr integer :: i, j, k, ce - real(pReal) :: phiDot, mobility phi_current = x_scal !-------------------------------------------------------------------------------------------------- @@ -272,8 +272,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(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)) + vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k)) enddo; enddo; enddo call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field @@ -281,10 +280,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call damage_nonlocal_getSourceAndItsTangent(phiDot, phi_current(i,j,k),ce) - mobility = damage_nonlocal_getMobility(ce) - scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) & - + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + mu_ref*phi_current(i,j,k) enddo; enddo; enddo @@ -293,6 +290,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) call utilities_FFTscalarForward call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc) call utilities_FFTscalarBackward + where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) & scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) & @@ -308,18 +306,18 @@ end subroutine formResidual !-------------------------------------------------------------------------------------------------- !> @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 mu_ref = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - ce = ce + 1 - K_ref = K_ref + damage_nonlocal_getDiffusion(ce) - mu_ref = mu_ref + damage_nonlocal_getMobility(ce) - enddo; enddo; enddo + do ce = 1, product(grid(1:2))*grid3 + K_ref = K_ref + homogenization_K_phi(ce) + mu_ref = mu_ref + homogenization_mu_phi(ce) + enddo + K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) mu_ref = mu_ref*wgt diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 14e2affb9..1ccb3b901 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -18,7 +18,6 @@ module grid_thermal_spectral use homogenization use YAML_types use config - use material implicit none private @@ -46,7 +45,7 @@ module grid_thermal_spectral !-------------------------------------------------------------------------------------------------- ! 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) :: mu_ref @@ -68,7 +67,7 @@ subroutine grid_thermal_spectral_init(T_0) PetscInt, dimension(0:worldsize-1) :: localK integer :: i, j, k, ce DM :: thermal_grid - PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscScalar, dimension(:,:,:), pointer :: T_PETSc PetscErrorCode :: ierr class(tNode), pointer :: & 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_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%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%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_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_thermal_rtol') !-------------------------------------------------------------------------------------------------- ! 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_lastInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal) + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(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) call homogenization_thermal_setField(T_0,0.0_pReal,ce) 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 DMDAVecRestoreArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) + + call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,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 @@ -190,9 +191,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call homogenization_thermal_setField(T_current(i,j,k), & - (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & - ce) + call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce) enddo; enddo; enddo call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) @@ -223,16 +222,14 @@ subroutine grid_thermal_spectral_forward(cutBack) !-------------------------------------------------------------------------------------------------- ! reverting thermal field state - ce = 0 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 x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current 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) ce = ce + 1 - call homogenization_thermal_setField(T_current(i,j,k), & - (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & - ce) + call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce) enddo; enddo; enddo else T_lastInc = T_current @@ -258,7 +255,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) PetscObject :: dummy PetscErrorCode :: ierr integer :: i, j, k, ce - real(pReal) :: Tdot T_current = x_scal !-------------------------------------------------------------------------------------------------- @@ -271,8 +267,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(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)) + vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k)) enddo; enddo; enddo call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field @@ -280,9 +275,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call thermal_conduction_getSource(Tdot,1,ce) - scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & - + homogenization_thermal_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & + + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + mu_ref*T_current(i,j,k) enddo; enddo; enddo @@ -302,18 +296,18 @@ end subroutine formResidual !-------------------------------------------------------------------------------------------------- !> @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 mu_ref = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - ce = ce + 1 - K_ref = K_ref + thermal_conduction_getConductivity(ce) - mu_ref = mu_ref + homogenization_thermal_mu_T(ce) - enddo; enddo; enddo + do ce = 1, product(grid(1:2))*grid3 + K_ref = K_ref + homogenization_K_T(ce) + mu_ref = mu_ref + homogenization_mu_T(ce) + enddo + K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) mu_ref = mu_ref*wgt diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 4ce75feca..3fadf1e96 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -100,10 +100,6 @@ module homogenization integer, intent(in) :: ce 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) real(pReal), intent(in) :: dt integer, intent(in) :: & @@ -135,46 +131,41 @@ module homogenization logical, dimension(2) :: doneAndHappy 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 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 - real(pReal) :: mu_T - end function homogenization_thermal_mu_T + real(pReal) :: f + end function homogenization_f_T module subroutine homogenization_thermal_setField(T,dot_T, ce) integer, intent(in) :: ce real(pReal), intent(in) :: T, dot_T 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 - real(pReal) :: T - end function homogenization_thermal_T + real(pReal) :: mu + end function homogenization_mu_phi - module subroutine thermal_conduction_getSource(Tdot, ip, el) - integer, intent(in) :: & - ip, & - el - real(pReal), intent(out) :: Tdot - end subroutine thermal_conduction_getSource - - module function damage_nonlocal_getMobility(ce) result(M) + module function homogenization_K_phi(ce) result(K) integer, intent(in) :: ce - real(pReal) :: M - end function damage_nonlocal_getMobility + real(pReal), dimension(3,3) :: K + 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 - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - phiDot - end subroutine damage_nonlocal_getSourceAndItsTangent + real(pReal), intent(in) :: phi + real(pReal) :: f + end function homogenization_f_phi module subroutine homogenization_set_phi(phi,ce) integer, intent(in) :: ce @@ -187,14 +178,14 @@ module homogenization public :: & homogenization_init, & materialpoint_stressAndItsTangent, & - homogenization_thermal_mu_T, & - thermal_conduction_getConductivity, & - thermal_conduction_getSource, & - damage_nonlocal_getMobility, & - damage_nonlocal_getSourceAndItsTangent, & - homogenization_set_phi, & + homogenization_mu_T, & + homogenization_K_T, & + homogenization_f_T, & homogenization_thermal_setfield, & - homogenization_thermal_T, & + homogenization_mu_phi, & + homogenization_K_phi, & + homogenization_f_phi, & + homogenization_set_phi, & homogenization_forward, & homogenization_results, & homogenization_restartRead, & @@ -202,9 +193,6 @@ module homogenization THERMAL_CONDUCTION_ID, & DAMAGE_NONLOCAL_ID - public :: & - damage_nonlocal_getDiffusion - contains @@ -314,7 +302,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE terminallyIll = .true. ! ...and kills all others endif enddo - call thermal_homogenize(ip,el) enddo enddo !$OMP END DO @@ -447,32 +434,6 @@ subroutine homogenization_restartRead(fileHandle) 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 !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index a067dde43..0546834fd 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -26,8 +26,8 @@ submodule(homogenization) damage type(tparameters), dimension(:), allocatable :: & param -contains +contains !-------------------------------------------------------------------------------------------------- !> @brief Allocate variables and set parameters. @@ -38,14 +38,12 @@ module subroutine damage_init() configHomogenizations, & configHomogenization, & configHomogenizationDamage, & - num_generic, & - material_homogenization - integer :: ho - integer :: Ninstances,Nmaterialpoints,h + num_generic + integer :: ho,Nmaterialpoints print'(/,a)', ' <<<+- homogenization:damage init -+>>>' - print'(/,a)', ' <<<+- homogenization:damage:pass init -+>>>' + configHomogenizations => config_material%get('homogenization') allocate(param(configHomogenizations%length)) @@ -62,6 +60,10 @@ module subroutine damage_init() #else prm%output = configHomogenizationDamage%get_as1dString('output',defaultVal=emptyStringArray) #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 prm%output = emptyStringArray endif @@ -73,18 +75,7 @@ module subroutine damage_init() num_generic => config_numerics%get('generic',defaultVal= emptyDict) num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) - Ninstances = count(damage_type == DAMAGE_nonlocal_ID) - - 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 + call pass_init() end subroutine damage_init @@ -100,53 +91,70 @@ module subroutine damage_partition(ce) if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return 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 - !-------------------------------------------------------------------------------------------------- -!> @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 - 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 real(pReal), intent(in) :: & phi - real(pReal), intent(out) :: & - phiDot + real(pReal) :: f - 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) integer, intent(in) :: ce real(pReal), intent(in) :: & phi + integer :: & ho, & en + ho = material_homogenizationID(ce) en = material_homogenizationEntry(ce) damagestate_h(ho)%state(1,en) = phi @@ -166,13 +174,13 @@ module subroutine damage_results(ho,group) integer :: o associate(prm => param(ho)) - outputsLoop: do o = 1,size(prm%output) - select case(prm%output(o)) - case ('phi') - call results_writeDataset(group,damagestate_h(ho)%state(1,:),prm%output(o),& - 'damage indicator','-') - end select - enddo outputsLoop + outputsLoop: do o = 1,size(prm%output) + select case(prm%output(o)) + case ('phi') + call results_writeDataset(group,damagestate_h(ho)%state(1,:),prm%output(o),& + 'damage indicator','-') + end select + enddo outputsLoop end associate end subroutine damage_results diff --git a/src/homogenization_damage_pass.f90 b/src/homogenization_damage_pass.f90 index cbb7ec79f..bdb44b822 100644 --- a/src/homogenization_damage_pass.f90 +++ b/src/homogenization_damage_pass.f90 @@ -6,8 +6,9 @@ submodule(homogenization:damage) damage_pass contains -module subroutine pass_init - +module subroutine pass_init() + + print'(/,a)', ' <<<+- homogenization:damage:pass init -+>>>' end subroutine pass_init diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index 7d1c64445..7ea1f74e9 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -32,25 +32,6 @@ submodule(homogenization) mechanical 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) logical, dimension(2) :: doneAndHappy real(pReal), dimension(:,:,:), intent(in) :: & @@ -148,39 +129,21 @@ module subroutine mechanical_homogenize(dt,ce) integer, intent(in) :: ce 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) = phase_P(1,ce) - homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) - - 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 + homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) & + / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) & + / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) end subroutine mechanical_homogenize diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 772d1c4f5..745c266d4 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -89,7 +89,8 @@ module subroutine RGC_init(num_homogMech) 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*, '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 associate(prm => param(material_homogenizationID(ce))) - + ho = material_homogenizationID(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 -!-------------------------------------------------------------------------------------------------- -!> @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 !-------------------------------------------------------------------------------------------------- @@ -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 = matmul(dst%orientation(1:3,1:3,en),interfaceNormal) ! map the normal vector into sample coordinate system (basis) - + end associate end function interfaceNormal diff --git a/src/homogenization_mechanical_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90 index 9a3704575..7b114d04f 100644 --- a/src/homogenization_mechanical_isostrain.f90 +++ b/src/homogenization_mechanical_isostrain.f90 @@ -6,21 +6,6 @@ !-------------------------------------------------------------------------------------------------- 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 !-------------------------------------------------------------------------------------------------- @@ -29,42 +14,21 @@ contains module subroutine isostrain_init integer :: & - h, & + ho, & Nmaterialpoints - class(tNode), pointer :: & - material_homogenization, & - homog, & - homogMech 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') - allocate(param(material_homogenization%length)) ! one container of parameters per homog + do ho = 1, size(homogenization_type) + if (homogenization_type(ho) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle - do h = 1, size(homogenization_type) - if (homogenization_type(h) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle - homog => material_homogenization%get(h) - homogMech => homog%get('mechanical') - 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 + Nmaterialpoints = count(material_homogenizationAt == ho) + homogState(ho)%sizeState = 0 + allocate(homogState(ho)%state0(0,Nmaterialpoints)) + allocate(homogState(ho)%state (0,Nmaterialpoints)) 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 + F = spread(avgF,3,size(F,3)) 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 diff --git a/src/homogenization_mechanical_pass.f90 b/src/homogenization_mechanical_pass.f90 index 23fe74f44..e2e44658a 100644 --- a/src/homogenization_mechanical_pass.f90 +++ b/src/homogenization_mechanical_pass.f90 @@ -14,25 +14,24 @@ contains module subroutine pass_init integer :: & - Ninstances, & - h, & + ho, & Nmaterialpoints print'(/,a)', ' <<<+- homogenization:mechanical:pass init -+>>>' - Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) + print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_NONE_ID) + flush(IO_STDOUT) - do h = 1, size(homogenization_type) - if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle + do ho = 1, size(homogenization_type) + 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)') - Nmaterialpoints = count(material_homogenizationAt == h) - homogState(h)%sizeState = 0 - allocate(homogState(h)%state0 (0,Nmaterialpoints)) - allocate(homogState(h)%state (0,Nmaterialpoints)) + Nmaterialpoints = count(material_homogenizationAt == ho) + homogState(ho)%sizeState = 0 + allocate(homogState(ho)%state0(0,Nmaterialpoints)) + allocate(homogState(ho)%state (0,Nmaterialpoints)) enddo diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 8dcddda7a..f681036de 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -45,8 +45,6 @@ module subroutine thermal_init() print'(/,a)', ' <<<+- homogenization:thermal init -+>>>' - print'(/,a)', ' <<<+- homogenization:thermal:pass init -+>>>' - configHomogenizations => config_material%get('homogenization') @@ -71,6 +69,8 @@ module subroutine thermal_init() end associate enddo + call pass_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 real(pReal), dimension(3,3) :: K - integer :: & - co + integer :: co - K = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) - K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseID(co,ce))) + K = phase_K_T(1,ce) + do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) + K = K + phase_K_T(co,ce) enddo K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) -end function thermal_conduction_getConductivity - - -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 +end function homogenization_K_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 - real(pReal) :: c_P + real(pReal) :: f 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)) - 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 - 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 -!-------------------------------------------------------------------------------------------------- -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) +!> @brief Set thermal field and its rate (T and dot_T). !-------------------------------------------------------------------------------------------------- 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 - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- @@ -219,44 +194,4 @@ module subroutine thermal_results(ho,group) 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 diff --git a/src/homogenization_thermal_pass.f90 b/src/homogenization_thermal_pass.f90 index 940a15024..87020d8d2 100644 --- a/src/homogenization_thermal_pass.f90 +++ b/src/homogenization_thermal_pass.f90 @@ -6,8 +6,9 @@ submodule(homogenization:thermal) thermal_pass contains -module subroutine pass_init - +module subroutine pass_init() + + print'(/,a)', ' <<<+- homogenization:thermal:pass init -+>>>' end subroutine pass_init diff --git a/src/lattice.f90 b/src/lattice.f90 index 5ec365fa7..06986b42b 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -394,13 +394,11 @@ module lattice ! SHOULD NOT BE PART OF LATTICE BEGIN real(pReal), dimension(:), allocatable, public, protected :: & lattice_mu, lattice_nu, & - lattice_M, & + lattice_mu_phi, & lattice_rho, & lattice_c_p real(pReal), dimension(:,:,:), allocatable, public, protected :: & - lattice_C66, & - lattice_K, & - lattice_D + lattice_C66 integer(kind(lattice_UNDEFINED_ID)), dimension(:), allocatable, public, protected :: & lattice_structure ! SHOULD NOT BE PART OF LATTICE END @@ -458,23 +456,19 @@ subroutine lattice_init phases, & phase, & mech, & - elasticity, & - thermal, & - damage + elasticity print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) +! SHOULD NOT BE PART OF LATTICE BEGIN + phases => config_material%get('phase') Nphases = phases%length allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID) allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) - allocate(lattice_K (3,3,Nphases), source=0.0_pReal) - allocate(lattice_D (3,3,Nphases), source=0.0_pReal) - - allocate(lattice_M,& - lattice_rho,lattice_c_p, & + allocate(lattice_rho, & lattice_mu, lattice_nu,& source=[(0.0_pReal,i=1,Nphases)]) @@ -522,32 +516,7 @@ subroutine lattice_init enddo lattice_rho(ph) = phase%get_asFloat('rho', defaultVal=0.0_pReal) - - ! 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 +! SHOULD NOT BE PART OF LATTICE END call selfTest diff --git a/src/phase.f90 b/src/phase.f90 index 1adba03d9..54116f3c9 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -145,26 +145,22 @@ module phase real(pReal), dimension(3,3) :: 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) integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: 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) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: 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) integer, intent(in) :: ph,me real(pReal) :: T @@ -188,13 +184,35 @@ module phase module subroutine phase_thermal_setField(T,dot_T, co,ce) real(pReal), intent(in) :: T, dot_T - integer, intent(in) :: ce, co + integer, intent(in) :: co, ce 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 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 =================================================================================== @@ -227,19 +245,18 @@ module phase 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 real(pReal), intent(in) :: & phi !< damage parameter real(pReal) :: & - phi_dot - end function phase_damage_phi_dot + f + 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 - real(pReal), intent(out) :: & - TDot - end subroutine phase_thermal_getRate + real(pReal) :: f + end function phase_f_T module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) integer, intent(in) :: & @@ -300,8 +317,12 @@ module phase public :: & phase_init, & phase_homogenizedC, & - phase_damage_phi_dot, & - phase_thermal_getRate, & + phase_f_phi, & + phase_f_T, & + phase_K_phi, & + phase_K_T, & + phase_mu_phi, & + phase_mu_T, & phase_results, & phase_allocateState, & phase_forward, & @@ -318,8 +339,7 @@ module phase phase_restartRead, & integrateDamageState, & phase_thermal_setField, & - phase_damage_set_phi, & - phase_damage_get_phi, & + phase_set_phi, & phase_P, & phase_set_F, & phase_F diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 62ce70368..8df2cd5e4 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -2,6 +2,12 @@ !> @brief internal microstructure state for all damage sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- 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 :: & DAMAGE_UNDEFINED_ID, & DAMAGE_ISOBRITTLE_ID, & @@ -16,13 +22,12 @@ submodule(phase) damage end type tDataContainer integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & - phase_source !< active sources mechanisms of each phase - - integer, dimension(:), allocatable :: & - phase_Nsources + phase_source !< active sources mechanisms of each phase type(tDataContainer), dimension(:), allocatable :: current + type(tDamageParameters), dimension(:), allocatable :: param + interface module function anisobrittle_init() result(mySources) @@ -100,18 +105,19 @@ module subroutine damage_init class(tNode), pointer :: & phases, & phase, & - sources - + sources, & + source + logical:: damage_active print'(/,a)', ' <<<+- phase:damage init -+>>>' phases => config_material%get('phase') allocate(current(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 Nmembers = count(material_phaseID == ph) @@ -122,14 +128,21 @@ module subroutine damage_init phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) 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 allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID) -! initialize source mechanisms - if(maxval(phase_Nsources) /= 0) then + if (damage_active) then where(isobrittle_init() ) phase_source = DAMAGE_ISOBRITTLE_ID where(isoductile_init() ) phase_source = DAMAGE_ISODUCTILE_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 !---------------------------------------------------------------------------------------------- -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 real(pReal), intent(in) :: & phi !< damage parameter real(pReal) :: & - phi_dot + f integer :: & ph, & @@ -159,13 +172,13 @@ module function phase_damage_phi_dot(phi,co,ce) result(phi_dot) select case(phase_source(ph)) case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID,DAMAGE_ANISODUCTILE_ID) - phi_dot = 1.0_pReal & - - phi*damageState(ph)%state(1,en) + f = 1.0_pReal & + - phi*damageState(ph)%state(1,en) case default - phi_dot = 0.0_pReal + f = 0.0_pReal 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 integer, intent(in) :: ph - integer :: so - - sourceLoop: do so = 1, phase_Nsources(ph) if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) & call results_closeGroup(results_addGroup(group//'damage')) - sourceType: select case (phase_source(ph)) + sourceType: select case (phase_source(ph)) - case (DAMAGE_ISOBRITTLE_ID) sourceType - call isobrittle_results(ph,group//'damage/') + case (DAMAGE_ISOBRITTLE_ID) sourceType + call isobrittle_results(ph,group//'damage/') - case (DAMAGE_ISODUCTILE_ID) sourceType - call isoductile_results(ph,group//'damage/') + case (DAMAGE_ISODUCTILE_ID) sourceType + call isoductile_results(ph,group//'damage/') - case (DAMAGE_ANISOBRITTLE_ID) sourceType - call anisobrittle_results(ph,group//'damage/') + case (DAMAGE_ANISOBRITTLE_ID) sourceType + call anisobrittle_results(ph,group//'damage/') - case (DAMAGE_ANISODUCTILE_ID) sourceType - call anisoductile_results(ph,group//'damage/') + case (DAMAGE_ANISODUCTILE_ID) sourceType + call anisoductile_results(ph,group//'damage/') - end select sourceType - - enddo SourceLoop + end select sourceType end subroutine damage_results @@ -339,6 +347,33 @@ function phase_damage_collectDotState(ph,me) result(broken) 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 @@ -362,19 +397,19 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) if (damageState(ph)%sizeState == 0) return - sourceType: select case (phase_source(ph)) + sourceType: select case (phase_source(ph)) - case (DAMAGE_ISOBRITTLE_ID) sourceType - call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) - broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me))) - if(.not. broken) then - myOffset = damageState(ph)%offsetDeltaState - mySize = damageState(ph)%sizeDeltaState - damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = & - damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me) - endif + case (DAMAGE_ISOBRITTLE_ID) sourceType + call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) + broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me))) + if(.not. broken) then + myOffset = damageState(ph)%offsetDeltaState + mySize = damageState(ph)%sizeDeltaState + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = & + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me) + endif - end select sourceType + end select sourceType end function phase_damage_deltaState @@ -411,7 +446,7 @@ end function source_active !---------------------------------------------------------------------------------------------- !< @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 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 -end subroutine phase_damage_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 +end subroutine phase_set_phi module function damage_phi(ph,me) result(phi) diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 1e157f338..f6d046266 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -3,6 +3,11 @@ !---------------------------------------------------------------------------------------------------- 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 :: & thermal_Nsources @@ -21,7 +26,9 @@ submodule(phase) thermal integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & 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 @@ -45,21 +52,19 @@ submodule(phase) thermal me end subroutine externalheat_dotState - module subroutine dissipation_getRate(TDot, ph,me) + module function dissipation_f_T(ph,me) result(f_T) integer, intent(in) :: & ph, & me - real(pReal), intent(out) :: & - TDot - end subroutine dissipation_getRate + real(pReal) :: f_T + end function dissipation_f_T - module subroutine externalheat_getRate(TDot, ph,me) + module function externalheat_f_T(ph,me) result(f_T) integer, intent(in) :: & ph, & me - real(pReal), intent(out) :: & - TDot - end subroutine externalheat_getRate + real(pReal) :: f_T + end function externalheat_f_T end interface @@ -87,6 +92,7 @@ module subroutine thermal_init(phases) allocate(thermalState(phases%length)) allocate(thermal_Nsources(phases%length),source = 0) + allocate(param(phases%length)) do ph = 1, phases%length Nmembers = count(material_phaseID == ph) @@ -94,9 +100,17 @@ module subroutine thermal_init(phases) allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal) phase => phases%get(ph) thermal => phase%get('thermal',defaultVal=emptyDict) + param(ph)%C_p = thermal%get_asFloat('c_p',defaultVal=0.0_pReal) + 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) thermal_Nsources(ph) = sources%length allocate(thermalstate(ph)%p(thermal_Nsources(ph))) + enddo 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 !---------------------------------------------------------------------------------------------- -module subroutine phase_thermal_getRate(TDot, ph,me) +module function phase_f_T(ph,me) result(f) integer, intent(in) :: ph, me - real(pReal), intent(out) :: & - TDot - - real(pReal) :: & - my_Tdot - integer :: & - so + real(pReal) :: f - TDot = 0.0_pReal + integer :: so + + + f = 0.0_pReal do so = 1, thermal_Nsources(ph) select case(thermal_source(so,ph)) + case (THERMAL_DISSIPATION_ID) - call dissipation_getRate(my_Tdot, ph,me) + f = f + dissipation_f_T(ph,me) 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 - Tdot = Tdot + my_Tdot + 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 +!-------------------------------------------------------------------------------------------------- +!> @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" ?? real(pReal), intent(in) :: Delta_t diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 index b16ed3cf7..3a4ee651a 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -69,17 +69,17 @@ end function dissipation_init !-------------------------------------------------------------------------------------------------- !> @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 - real(pReal), intent(out) :: & - TDot + real(pReal) :: & + f_T 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 subroutine dissipation_getRate +end function dissipation_f_T end submodule dissipation diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index d5bbc7c38..6d4403ab8 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -100,13 +100,13 @@ end subroutine externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @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) :: & ph, & me - real(pReal), intent(out) :: & - TDot + real(pReal) :: & + f_T integer :: & so, interval @@ -122,12 +122,12 @@ module subroutine externalheat_getRate(TDot, ph, me) if ( (frac_time < 0.0_pReal .and. interval == 1) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & - TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + & - prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... + f_T = prm%f_T(interval ) * (1.0_pReal - frac_time) + & + prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... ! ...or extrapolate if outside me bounds enddo end associate -end subroutine externalheat_getRate +end function externalheat_f_T end submodule externalheat