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

Merge branch 'development' into 29-rename-mesh_element-array

This commit is contained in:
Martin Diehl 2018-10-01 22:29:31 +02:00
parent 49a4202e26 27bde05529
commit 98296a978e
31 changed files with 325378 additions and 1302 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
.gitlab-ci.yml
View File

@ -60,18 +60,17 @@ variables:
IntelCompiler: "$IntelCompiler18_1" IntelCompiler: "$IntelCompiler18_1"
GNUCompiler: "$GNUCompiler7_3" GNUCompiler: "$GNUCompiler7_3"
# ++++++++++++ MPI +++++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++ MPI +++++++++++++++++++++++++++++++++++++++++++++++++++
MPICH3_2Intel17_0: "MPI/Intel/17.0/MPICH/3.2"
MPICH3_2Intel18_1: "MPI/Intel/18.1/MPICH/3.2.1" MPICH3_2Intel18_1: "MPI/Intel/18.1/MPICH/3.2.1"
MPICH3_2GNU7_3: "MPI/GNU/7.3/MPICH/3.2.1" MPICH3_2GNU7_3: "MPI/GNU/7.3/MPICH/3.2.1"
# ------------ Defaults ---------------------------------------------- # ------------ Defaults ----------------------------------------------
MPICH_GNU: "$MPICH3_2GNU7_3"
MPICH_Intel: "$MPICH3_2Intel18_1" MPICH_Intel: "$MPICH3_2Intel18_1"
MPICH_GNU: "$MPICH3_2GNU7_3"
# ++++++++++++ PETSc +++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++ PETSc +++++++++++++++++++++++++++++++++++++++++++++++++
PETSc3_9_1MPICH3_2Intel18_1: "Libraries/PETSc/3.9.1/Intel-18.1-MPICH-3.2.1" PETSc3_10_0MPICH3_2Intel18_1: "Libraries/PETSc/3.10.0/Intel-18.1-MPICH-3.2.1"
PETSc3_9_1MPICH3_2GNU7_3: "Libraries/PETSc/3.9.1/GNU-7.3-MPICH-3.2.1" PETSc3_10_0MPICH3_2GNU7_3: "Libraries/PETSc/3.10.0/GNU-7.3-MPICH-3.2.1"
# ------------ Defaults ---------------------------------------------- # ------------ Defaults ----------------------------------------------
PETSc_MPICH_Intel: "$PETSc3_9_1MPICH3_2Intel18_1" PETSc_MPICH_Intel: "$PETSc3_10_0MPICH3_2Intel18_1"
PETSc_MPICH_GNU: "$PETSc3_9_1MPICH3_2GNU7_3" PETSc_MPICH_GNU: "$PETSc3_10_0MPICH3_2GNU7_3"
# ++++++++++++ FEM +++++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++ FEM +++++++++++++++++++++++++++++++++++++++++++++++++++
Abaqus2017: "FEM/Abaqus/2017" Abaqus2017: "FEM/Abaqus/2017"
MSC2017: "FEM/MSC/2017" MSC2017: "FEM/MSC/2017"

2
PRIVATE

@ -1 +1 @@
Subproject commit aeae4513b1ffb43b949399c12bae27fc6abb6f29 Subproject commit b893d0488ec211d4d2cb02c284079e6fbcd89baa

2
VERSION
View File

@ -1 +1 @@
v2.0.2-582-g0710609c v2.0.2-591-ga00d15b8

8
examples/.gitignore vendored
View File

@ -1 +1,9 @@
*.C_ref
*.mesh
*.outputConstitutive
*.outputCrystallite
*.outputHomogenization
*.spectralOut
*.sta
*.vt*
postProc postProc

759
examples/FEM/polyXtal/material.config Normal file
View File

@ -0,0 +1,759 @@
#-------------------#
<homogenization>
#-------------------#
[SX]
mech none
#-------------------#
<crystallite>
#-------------------#
[aLittleSomething]
(output) f
(output) p
#-------------------#
<phase>
#-------------------#
[Aluminum]
elasticity hooke
plasticity phenopowerlaw
lattice_structure fcc
Nslip 12 # per family
Ntwin 0 # per family
c11 106.75e9
c12 60.41e9
c44 28.34e9
gdot0_slip 0.001
n_slip 20
tau0_slip 31e6 # per family
tausat_slip 63e6 # per family
a_slip 2.25
gdot0_twin 0.001
n_twin 20
tau0_twin 31e6 # per family
s_pr 0 # push-up factor for slip saturation due to twinning
twin_b 0
twin_c 0
twin_d 0
twin_e 0
h0_slipslip 75e6
h0_twinslip 0
h0_twintwin 0
interaction_slipslip 1 1 1.4 1.4 1.4 1.4
interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
atol_resistance 1
#-------------------#
<microstructure>
#-------------------#
[Grain001]
crystallite 1
(constituent) phase 1 texture 1 fraction 1.0
[Grain002]
crystallite 1
(constituent) phase 1 texture 2 fraction 1.0
[Grain003]
crystallite 1
(constituent) phase 1 texture 3 fraction 1.0
[Grain004]
crystallite 1
(constituent) phase 1 texture 4 fraction 1.0
[Grain005]
crystallite 1
(constituent) phase 1 texture 5 fraction 1.0
[Grain006]
crystallite 1
(constituent) phase 1 texture 6 fraction 1.0
[Grain007]
crystallite 1
(constituent) phase 1 texture 7 fraction 1.0
[Grain008]
crystallite 1
(constituent) phase 1 texture 8 fraction 1.0
[Grain009]
crystallite 1
(constituent) phase 1 texture 9 fraction 1.0
[Grain010]
crystallite 1
(constituent) phase 1 texture 10 fraction 1.0
[Grain011]
crystallite 1
(constituent) phase 1 texture 11 fraction 1.0
[Grain012]
crystallite 1
(constituent) phase 1 texture 12 fraction 1.0
[Grain013]
crystallite 1
(constituent) phase 1 texture 13 fraction 1.0
[Grain014]
crystallite 1
(constituent) phase 1 texture 14 fraction 1.0
[Grain015]
crystallite 1
(constituent) phase 1 texture 15 fraction 1.0
[Grain016]
crystallite 1
(constituent) phase 1 texture 16 fraction 1.0
[Grain017]
crystallite 1
(constituent) phase 1 texture 17 fraction 1.0
[Grain018]
crystallite 1
(constituent) phase 1 texture 18 fraction 1.0
[Grain019]
crystallite 1
(constituent) phase 1 texture 19 fraction 1.0
[Grain020]
crystallite 1
(constituent) phase 1 texture 20 fraction 1.0
[Grain021]
crystallite 1
(constituent) phase 1 texture 21 fraction 1.0
[Grain022]
crystallite 1
(constituent) phase 1 texture 22 fraction 1.0
[Grain023]
crystallite 1
(constituent) phase 1 texture 23 fraction 1.0
[Grain024]
crystallite 1
(constituent) phase 1 texture 24 fraction 1.0
[Grain025]
crystallite 1
(constituent) phase 1 texture 25 fraction 1.0
[Grain026]
crystallite 1
(constituent) phase 1 texture 26 fraction 1.0
[Grain027]
crystallite 1
(constituent) phase 1 texture 27 fraction 1.0
[Grain028]
crystallite 1
(constituent) phase 1 texture 28 fraction 1.0
[Grain029]
crystallite 1
(constituent) phase 1 texture 29 fraction 1.0
[Grain030]
crystallite 1
(constituent) phase 1 texture 30 fraction 1.0
[Grain031]
crystallite 1
(constituent) phase 1 texture 31 fraction 1.0
[Grain032]
crystallite 1
(constituent) phase 1 texture 32 fraction 1.0
[Grain033]
crystallite 1
(constituent) phase 1 texture 33 fraction 1.0
[Grain034]
crystallite 1
(constituent) phase 1 texture 34 fraction 1.0
[Grain035]
crystallite 1
(constituent) phase 1 texture 35 fraction 1.0
[Grain036]
crystallite 1
(constituent) phase 1 texture 36 fraction 1.0
[Grain037]
crystallite 1
(constituent) phase 1 texture 37 fraction 1.0
[Grain038]
crystallite 1
(constituent) phase 1 texture 38 fraction 1.0
[Grain039]
crystallite 1
(constituent) phase 1 texture 39 fraction 1.0
[Grain040]
crystallite 1
(constituent) phase 1 texture 40 fraction 1.0
[Grain041]
crystallite 1
(constituent) phase 1 texture 41 fraction 1.0
[Grain042]
crystallite 1
(constituent) phase 1 texture 42 fraction 1.0
[Grain043]
crystallite 1
(constituent) phase 1 texture 43 fraction 1.0
[Grain044]
crystallite 1
(constituent) phase 1 texture 44 fraction 1.0
[Grain045]
crystallite 1
(constituent) phase 1 texture 45 fraction 1.0
[Grain046]
crystallite 1
(constituent) phase 1 texture 46 fraction 1.0
[Grain047]
crystallite 1
(constituent) phase 1 texture 47 fraction 1.0
[Grain048]
crystallite 1
(constituent) phase 1 texture 48 fraction 1.0
[Grain049]
crystallite 1
(constituent) phase 1 texture 49 fraction 1.0
[Grain050]
crystallite 1
(constituent) phase 1 texture 50 fraction 1.0
[Grain051]
crystallite 1
(constituent) phase 1 texture 51 fraction 1.0
[Grain052]
crystallite 1
(constituent) phase 1 texture 52 fraction 1.0
[Grain053]
crystallite 1
(constituent) phase 1 texture 53 fraction 1.0
[Grain054]
crystallite 1
(constituent) phase 1 texture 54 fraction 1.0
[Grain055]
crystallite 1
(constituent) phase 1 texture 55 fraction 1.0
[Grain056]
crystallite 1
(constituent) phase 1 texture 56 fraction 1.0
[Grain057]
crystallite 1
(constituent) phase 1 texture 57 fraction 1.0
[Grain058]
crystallite 1
(constituent) phase 1 texture 58 fraction 1.0
[Grain059]
crystallite 1
(constituent) phase 1 texture 59 fraction 1.0
[Grain060]
crystallite 1
(constituent) phase 1 texture 60 fraction 1.0
[Grain061]
crystallite 1
(constituent) phase 1 texture 61 fraction 1.0
[Grain062]
crystallite 1
(constituent) phase 1 texture 62 fraction 1.0
[Grain063]
crystallite 1
(constituent) phase 1 texture 63 fraction 1.0
[Grain064]
crystallite 1
(constituent) phase 1 texture 64 fraction 1.0
[Grain065]
crystallite 1
(constituent) phase 1 texture 65 fraction 1.0
[Grain066]
crystallite 1
(constituent) phase 1 texture 66 fraction 1.0
[Grain067]
crystallite 1
(constituent) phase 1 texture 67 fraction 1.0
[Grain068]
crystallite 1
(constituent) phase 1 texture 68 fraction 1.0
[Grain069]
crystallite 1
(constituent) phase 1 texture 69 fraction 1.0
[Grain070]
crystallite 1
(constituent) phase 1 texture 70 fraction 1.0
[Grain071]
crystallite 1
(constituent) phase 1 texture 71 fraction 1.0
[Grain072]
crystallite 1
(constituent) phase 1 texture 72 fraction 1.0
[Grain073]
crystallite 1
(constituent) phase 1 texture 73 fraction 1.0
[Grain074]
crystallite 1
(constituent) phase 1 texture 74 fraction 1.0
[Grain075]
crystallite 1
(constituent) phase 1 texture 75 fraction 1.0
[Grain076]
crystallite 1
(constituent) phase 1 texture 76 fraction 1.0
[Grain077]
crystallite 1
(constituent) phase 1 texture 77 fraction 1.0
[Grain078]
crystallite 1
(constituent) phase 1 texture 78 fraction 1.0
[Grain079]
crystallite 1
(constituent) phase 1 texture 79 fraction 1.0
[Grain080]
crystallite 1
(constituent) phase 1 texture 80 fraction 1.0
[Grain081]
crystallite 1
(constituent) phase 1 texture 81 fraction 1.0
[Grain082]
crystallite 1
(constituent) phase 1 texture 82 fraction 1.0
[Grain083]
crystallite 1
(constituent) phase 1 texture 83 fraction 1.0
[Grain084]
crystallite 1
(constituent) phase 1 texture 84 fraction 1.0
[Grain085]
crystallite 1
(constituent) phase 1 texture 85 fraction 1.0
[Grain086]
crystallite 1
(constituent) phase 1 texture 86 fraction 1.0
[Grain087]
crystallite 1
(constituent) phase 1 texture 87 fraction 1.0
[Grain088]
crystallite 1
(constituent) phase 1 texture 88 fraction 1.0
[Grain089]
crystallite 1
(constituent) phase 1 texture 89 fraction 1.0
[Grain090]
crystallite 1
(constituent) phase 1 texture 90 fraction 1.0
[Grain091]
crystallite 1
(constituent) phase 1 texture 91 fraction 1.0
[Grain092]
crystallite 1
(constituent) phase 1 texture 92 fraction 1.0
[Grain093]
crystallite 1
(constituent) phase 1 texture 93 fraction 1.0
[Grain094]
crystallite 1
(constituent) phase 1 texture 94 fraction 1.0
[Grain095]
crystallite 1
(constituent) phase 1 texture 95 fraction 1.0
[Grain096]
crystallite 1
(constituent) phase 1 texture 96 fraction 1.0
[Grain097]
crystallite 1
(constituent) phase 1 texture 97 fraction 1.0
[Grain098]
crystallite 1
(constituent) phase 1 texture 98 fraction 1.0
[Grain099]
crystallite 1
(constituent) phase 1 texture 99 fraction 1.0
[Grain100]
crystallite 1
(constituent) phase 1 texture 100 fraction 1.0
#-------------------#
<texture>
#-------------------#
[Grain001]
(gauss) phi1 172.344 Phi 114.046 phi2 294.669 scatter 0.0 fraction 1.0
[Grain002]
(gauss) phi1 186.013 Phi 94.7338 phi2 329.683 scatter 0.0 fraction 1.0
[Grain003]
(gauss) phi1 162.41 Phi 98.9455 phi2 130.322 scatter 0.0 fraction 1.0
[Grain004]
(gauss) phi1 355.272 Phi 140.621 phi2 125.567 scatter 0.0 fraction 1.0
[Grain005]
(gauss) phi1 21.7641 Phi 143.388 phi2 240.373 scatter 0.0 fraction 1.0
[Grain006]
(gauss) phi1 88.1966 Phi 92.3358 phi2 194.78 scatter 0.0 fraction 1.0
[Grain007]
(gauss) phi1 161.137 Phi 78.0062 phi2 111.948 scatter 0.0 fraction 1.0
[Grain008]
(gauss) phi1 169.792 Phi 89.5333 phi2 159.265 scatter 0.0 fraction 1.0
[Grain009]
(gauss) phi1 264.847 Phi 130.291 phi2 180.604 scatter 0.0 fraction 1.0
[Grain010]
(gauss) phi1 70.6323 Phi 84.1754 phi2 341.162 scatter 0.0 fraction 1.0
[Grain011]
(gauss) phi1 67.7751 Phi 36.1662 phi2 139.898 scatter 0.0 fraction 1.0
[Grain012]
(gauss) phi1 111.621 Phi 19.1089 phi2 228.338 scatter 0.0 fraction 1.0
[Grain013]
(gauss) phi1 129.9 Phi 139.011 phi2 238.735 scatter 0.0 fraction 1.0
[Grain014]
(gauss) phi1 221.405 Phi 129.743 phi2 99.6471 scatter 0.0 fraction 1.0
[Grain015]
(gauss) phi1 241.783 Phi 98.3729 phi2 260.615 scatter 0.0 fraction 1.0
[Grain016]
(gauss) phi1 72.5592 Phi 122.403 phi2 165.046 scatter 0.0 fraction 1.0
[Grain017]
(gauss) phi1 64.8818 Phi 82.6384 phi2 236.305 scatter 0.0 fraction 1.0
[Grain018]
(gauss) phi1 201.096 Phi 65.9312 phi2 330.745 scatter 0.0 fraction 1.0
[Grain019]
(gauss) phi1 192.994 Phi 81.9371 phi2 239.326 scatter 0.0 fraction 1.0
[Grain020]
(gauss) phi1 125.335 Phi 90.4527 phi2 207.982 scatter 0.0 fraction 1.0
[Grain021]
(gauss) phi1 55.8848 Phi 26.4455 phi2 100.921 scatter 0.0 fraction 1.0
[Grain022]
(gauss) phi1 40.722 Phi 95.6415 phi2 269.174 scatter 0.0 fraction 1.0
[Grain023]
(gauss) phi1 250.487 Phi 69.6035 phi2 201.732 scatter 0.0 fraction 1.0
[Grain024]
(gauss) phi1 204.199 Phi 84.983 phi2 20.3469 scatter 0.0 fraction 1.0
[Grain025]
(gauss) phi1 12.7416 Phi 128.589 phi2 271.553 scatter 0.0 fraction 1.0
[Grain026]
(gauss) phi1 299.704 Phi 85.3961 phi2 217.359 scatter 0.0 fraction 1.0
[Grain027]
(gauss) phi1 48.8232 Phi 83.6209 phi2 200.361 scatter 0.0 fraction 1.0
[Grain028]
(gauss) phi1 336.395 Phi 97.3059 phi2 187.071 scatter 0.0 fraction 1.0
[Grain029]
(gauss) phi1 274.354 Phi 78.2424 phi2 320.308 scatter 0.0 fraction 1.0
[Grain030]
(gauss) phi1 320.776 Phi 149.72 phi2 163.862 scatter 0.0 fraction 1.0
[Grain031]
(gauss) phi1 179.549 Phi 106.548 phi2 345.498 scatter 0.0 fraction 1.0
[Grain032]
(gauss) phi1 163.508 Phi 24.4238 phi2 127.809 scatter 0.0 fraction 1.0
[Grain033]
(gauss) phi1 193.405 Phi 157.012 phi2 321.342 scatter 0.0 fraction 1.0
[Grain034]
(gauss) phi1 9.09886 Phi 95.9453 phi2 102.32 scatter 0.0 fraction 1.0
[Grain035]
(gauss) phi1 353.876 Phi 150.824 phi2 174.886 scatter 0.0 fraction 1.0
[Grain036]
(gauss) phi1 138.914 Phi 76.5811 phi2 167.927 scatter 0.0 fraction 1.0
[Grain037]
(gauss) phi1 262.655 Phi 76.2738 phi2 12.4459 scatter 0.0 fraction 1.0
[Grain038]
(gauss) phi1 121.849 Phi 65.5254 phi2 192.601 scatter 0.0 fraction 1.0
[Grain039]
(gauss) phi1 275.824 Phi 81.6788 phi2 164.228 scatter 0.0 fraction 1.0
[Grain040]
(gauss) phi1 68.9202 Phi 160.5 phi2 210.862 scatter 0.0 fraction 1.0
[Grain041]
(gauss) phi1 51.0398 Phi 82.7291 phi2 74.016 scatter 0.0 fraction 1.0
[Grain042]
(gauss) phi1 338.746 Phi 62.7854 phi2 129.362 scatter 0.0 fraction 1.0
[Grain043]
(gauss) phi1 204.51 Phi 151.256 phi2 178.89 scatter 0.0 fraction 1.0
[Grain044]
(gauss) phi1 122.098 Phi 104.003 phi2 323.04 scatter 0.0 fraction 1.0
[Grain045]
(gauss) phi1 106.693 Phi 108.61 phi2 336.935 scatter 0.0 fraction 1.0
[Grain046]
(gauss) phi1 118.856 Phi 160.992 phi2 316.152 scatter 0.0 fraction 1.0
[Grain047]
(gauss) phi1 177.962 Phi 114.868 phi2 13.6918 scatter 0.0 fraction 1.0
[Grain048]
(gauss) phi1 330.273 Phi 174.495 phi2 231.249 scatter 0.0 fraction 1.0
[Grain049]
(gauss) phi1 7.31937 Phi 94.7313 phi2 17.8461 scatter 0.0 fraction 1.0
[Grain050]
(gauss) phi1 74.3385 Phi 49.9546 phi2 286.482 scatter 0.0 fraction 1.0
[Grain051]
(gauss) phi1 326.388 Phi 76.9547 phi2 214.758 scatter 0.0 fraction 1.0
[Grain052]
(gauss) phi1 276.024 Phi 72.1242 phi2 275.884 scatter 0.0 fraction 1.0
[Grain053]
(gauss) phi1 137.681 Phi 116.99 phi2 6.87047 scatter 0.0 fraction 1.0
[Grain054]
(gauss) phi1 200.213 Phi 123.618 phi2 268.84 scatter 0.0 fraction 1.0
[Grain055]
(gauss) phi1 7.13702 Phi 56.2015 phi2 119.65 scatter 0.0 fraction 1.0
[Grain056]
(gauss) phi1 72.1783 Phi 81.0906 phi2 6.06213 scatter 0.0 fraction 1.0
[Grain057]
(gauss) phi1 184.565 Phi 110.01 phi2 239.546 scatter 0.0 fraction 1.0
[Grain058]
(gauss) phi1 210.124 Phi 128.631 phi2 8.61611 scatter 0.0 fraction 1.0
[Grain059]
(gauss) phi1 290.326 Phi 170.412 phi2 144.269 scatter 0.0 fraction 1.0
[Grain060]
(gauss) phi1 204.748 Phi 76.7343 phi2 200.385 scatter 0.0 fraction 1.0
[Grain061]
(gauss) phi1 54.3015 Phi 65.9143 phi2 117.373 scatter 0.0 fraction 1.0
[Grain062]
(gauss) phi1 261.263 Phi 52.255 phi2 95.9146 scatter 0.0 fraction 1.0
[Grain063]
(gauss) phi1 328.054 Phi 51.0778 phi2 24.2782 scatter 0.0 fraction 1.0
[Grain064]
(gauss) phi1 163.03 Phi 154.894 phi2 64.126 scatter 0.0 fraction 1.0
[Grain065]
(gauss) phi1 183.87 Phi 80.1848 phi2 18.7438 scatter 0.0 fraction 1.0
[Grain066]
(gauss) phi1 219.91 Phi 113.727 phi2 126.67 scatter 0.0 fraction 1.0
[Grain067]
(gauss) phi1 1.43844 Phi 87.6365 phi2 217.342 scatter 0.0 fraction 1.0
[Grain068]
(gauss) phi1 16.6245 Phi 162.07 phi2 43.7899 scatter 0.0 fraction 1.0
[Grain069]
(gauss) phi1 16.86 Phi 53.8682 phi2 256.917 scatter 0.0 fraction 1.0
[Grain070]
(gauss) phi1 1.01921 Phi 118.449 phi2 307.223 scatter 0.0 fraction 1.0
[Grain071]
(gauss) phi1 19.0397 Phi 83.8885 phi2 262.687 scatter 0.0 fraction 1.0
[Grain072]
(gauss) phi1 99.799 Phi 77.2307 phi2 84.9727 scatter 0.0 fraction 1.0
[Grain073]
(gauss) phi1 234.292 Phi 63.5029 phi2 250.315 scatter 0.0 fraction 1.0
[Grain074]
(gauss) phi1 315.529 Phi 106.015 phi2 103.711 scatter 0.0 fraction 1.0
[Grain075]
(gauss) phi1 235.595 Phi 110.152 phi2 210.277 scatter 0.0 fraction 1.0
[Grain076]
(gauss) phi1 341.907 Phi 17.1839 phi2 332.75 scatter 0.0 fraction 1.0
[Grain077]
(gauss) phi1 352.166 Phi 88.6049 phi2 114.964 scatter 0.0 fraction 1.0
[Grain078]
(gauss) phi1 342.33 Phi 117.777 phi2 180.346 scatter 0.0 fraction 1.0
[Grain079]
(gauss) phi1 224.952 Phi 70.5702 phi2 148.486 scatter 0.0 fraction 1.0
[Grain080]
(gauss) phi1 7.71702 Phi 23.6124 phi2 131.591 scatter 0.0 fraction 1.0
[Grain081]
(gauss) phi1 65.1024 Phi 138.774 phi2 247.344 scatter 0.0 fraction 1.0
[Grain082]
(gauss) phi1 37.6181 Phi 51.5209 phi2 8.4169 scatter 0.0 fraction 1.0
[Grain083]
(gauss) phi1 245.335 Phi 53.4543 phi2 52.5205 scatter 0.0 fraction 1.0
[Grain084]
(gauss) phi1 259.572 Phi 87.7026 phi2 272.065 scatter 0.0 fraction 1.0
[Grain085]
(gauss) phi1 269.39 Phi 103.379 phi2 132.506 scatter 0.0 fraction 1.0
[Grain086]
(gauss) phi1 175.156 Phi 119.338 phi2 355.51 scatter 0.0 fraction 1.0
[Grain087]
(gauss) phi1 248.11 Phi 39.4772 phi2 310.371 scatter 0.0 fraction 1.0
[Grain088]
(gauss) phi1 121.809 Phi 141.465 phi2 10.0736 scatter 0.0 fraction 1.0
[Grain089]
(gauss) phi1 2.4357 Phi 47.118 phi2 274.654 scatter 0.0 fraction 1.0
[Grain090]
(gauss) phi1 314.188 Phi 134.146 phi2 250.673 scatter 0.0 fraction 1.0
[Grain091]
(gauss) phi1 114.815 Phi 121.132 phi2 275.124 scatter 0.0 fraction 1.0
[Grain092]
(gauss) phi1 126.699 Phi 99.0325 phi2 320.537 scatter 0.0 fraction 1.0
[Grain093]
(gauss) phi1 184.138 Phi 20.1663 phi2 159.314 scatter 0.0 fraction 1.0
[Grain094]
(gauss) phi1 296.502 Phi 15.2389 phi2 39.382 scatter 0.0 fraction 1.0
[Grain095]
(gauss) phi1 167.8 Phi 151.764 phi2 192.568 scatter 0.0 fraction 1.0
[Grain096]
(gauss) phi1 257.822 Phi 133.446 phi2 257.108 scatter 0.0 fraction 1.0
[Grain097]
(gauss) phi1 71.6923 Phi 74.5726 phi2 342.575 scatter 0.0 fraction 1.0
[Grain098]
(gauss) phi1 176.748 Phi 28.39 phi2 327.375 scatter 0.0 fraction 1.0
[Grain099]
(gauss) phi1 121.822 Phi 141.836 phi2 22.6349 scatter 0.0 fraction 1.0
[Grain100]
(gauss) phi1 180.151 Phi 109.246 phi2 146.177 scatter 0.0 fraction 1.0

3
examples/FEM/polyXtal/numerics.config Normal file
View File

@ -0,0 +1,3 @@
residualStiffness 0.001
charLength 0.02
petsc_options -mech_snes_type newtonls -mech_ksp_type fgmres -mech_pc_type ml -mech_ksp_monitor

210313
examples/FEM/polyXtal/polyXtal_100grains.msh Normal file
View File

File diff suppressed because it is too large Load Diff

113546
examples/FEM/polyXtal/polyXtal_20grains.msh Normal file
View File

File diff suppressed because it is too large Load Diff

9
examples/FEM/polyXtal/tensionZ.load Normal file
View File

@ -0,0 +1,9 @@
#initial elastic step
$Loadcase 1 time 0.0005 incs 1 frequency 5
Face 1 X 0.01
Face 2 X 0.00
$EndLoadcase
$Loadcase 2 time 10.0 incs 200 frequency 5
Face 1 X 0.01
Face 2 X 0.00
$EndLoadcase

8
examples/SpectralMethod/.gitignore vendored
View File

@ -1,8 +0,0 @@
*.C_ref
*.mesh
*.outputConstitutive
*.outputCrystallite
*.outputHomogenization
*.spectralOut
*.sta
*.vtk

1
examples/SpectralMethod/Polycrystal/material.config
View File

@ -54,7 +54,6 @@ h0_slipslip 75e6
interaction_slipslip 1 1 1.4 1.4 1.4 1.4 interaction_slipslip 1 1 1.4 1.4 1.4 1.4
atol_resistance 1 atol_resistance 1
#-------------------# #-------------------#
<microstructure> <microstructure>
#-------------------# #-------------------#

79
lib/damask/asciitable.py
View File

@ -1,14 +1,14 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys import os, sys
import numpy as np import numpy as np
# ------------------------------------------------------------------ # ------------------------------------------------------------------
# python 3 has no unicode object, this ensures that the code works on Python 2&3 # python 3 has no unicode object, this ensures that the code works on Python 2&3
try: try:
test=isinstance('test', unicode) test = isinstance('test', unicode)
except(NameError): except NameError:
unicode=str unicode = str
# ------------------------------------------------------------------ # ------------------------------------------------------------------
class ASCIItable(): class ASCIItable():
@ -63,16 +63,17 @@ class ASCIItable():
x): x):
try: try:
return float(x) return float(x)
except: except ValueError:
return 0.0 return 0.0
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def _removeCRLF(self, def _removeCRLF(self,
string): string):
"""Delete any carriage return and line feed from string"""
try: try:
return string.replace('\n','').replace('\r','') return string.replace('\n','').replace('\r','')
except: except AttributeError:
return string return str(string)
# ------------------------------------------------------------------ # ------------------------------------------------------------------
@ -93,22 +94,22 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def input_close(self): def input_close(self):
try: # try:
if self.__IO__['in'] != sys.stdin: self.__IO__['in'].close() if self.__IO__['in'] != sys.stdin: self.__IO__['in'].close()
except: # except:
pass # pass
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def output_write(self, def output_write(self,
what): what):
"""Aggregate a single row (string) or list of (possibly containing further lists of) rows into output""" """Aggregate a single row (string) or list of (possibly containing further lists of) rows into output"""
if not isinstance(what, (str, unicode)): if isinstance(what, (str, unicode)):
self.__IO__['output'] += [what]
else:
try: try:
for item in what: self.output_write(item) for item in what: self.output_write(item)
except: except TypeError:
self.__IO__['output'] += [str(what)] self.__IO__['output'] += [str(what)]
else:
self.__IO__['output'] += [what]
return self.__IO__['buffered'] or self.output_flush() return self.__IO__['buffered'] or self.output_flush()
@ -129,10 +130,10 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def output_close(self, def output_close(self,
dismiss = False): dismiss = False):
try: # try:
if self.__IO__['out'] != sys.stdout: self.__IO__['out'].close() if self.__IO__['out'] != sys.stdout: self.__IO__['out'].close()
except: # except:
pass # pass
if dismiss and os.path.isfile(self.__IO__['out'].name): if dismiss and os.path.isfile(self.__IO__['out'].name):
os.remove(self.__IO__['out'].name) os.remove(self.__IO__['out'].name)
elif self.__IO__['inPlace']: elif self.__IO__['inPlace']:
@ -150,7 +151,7 @@ class ASCIItable():
try: try:
self.__IO__['in'].seek(0) self.__IO__['in'].seek(0)
except: except IOError:
pass pass
firstline = self.__IO__['in'].readline().strip() firstline = self.__IO__['in'].readline().strip()
@ -170,7 +171,7 @@ class ASCIItable():
else: # other table format else: # other table format
try: try:
self.__IO__['in'].seek(0) # try to rewind self.__IO__['in'].seek(0) # try to rewind
except: except IOError:
self.__IO__['readBuffer'] = [firstline] # or at least save data in buffer self.__IO__['readBuffer'] = [firstline] # or at least save data in buffer
while self.data_read(advance = False, respectLabels = False): while self.data_read(advance = False, respectLabels = False):
@ -197,7 +198,9 @@ class ASCIItable():
"""Write current header information (info + labels)""" """Write current header information (info + labels)"""
head = ['{}\theader'.format(len(self.info)+self.__IO__['labeled'])] if header else [] head = ['{}\theader'.format(len(self.info)+self.__IO__['labeled'])] if header else []
head.append(self.info) head.append(self.info)
if self.__IO__['labeled']: head.append('\t'.join(map(self._quote,self.tags))) if self.__IO__['labeled']:
head.append('\t'.join(map(self._quote,self.tags)))
if len(self.tags) == 0: raise ValueError('no labels present.')
return self.output_write(head) return self.output_write(head)
@ -257,13 +260,13 @@ class ASCIItable():
what, what,
reset = False): reset = False):
"""Add item or list to existing set of labels (and switch on labeling)""" """Add item or list to existing set of labels (and switch on labeling)"""
if not isinstance(what, (str, unicode)): if isinstance(what, (str, unicode)):
self.tags += [self._removeCRLF(what)]
else:
try: try:
for item in what: self.labels_append(item) for item in what: self.labels_append(item)
except: except TypeError:
self.tags += [self._removeCRLF(str(what))] self.tags += [self._removeCRLF(str(what))]
else:
self.tags += [self._removeCRLF(what)]
self.__IO__['labeled'] = True # switch on processing (in particular writing) of tags self.__IO__['labeled'] = True # switch on processing (in particular writing) of tags
if reset: self.__IO__['tags'] = list(self.tags) # subsequent data_read uses current tags as data size if reset: self.__IO__['tags'] = list(self.tags) # subsequent data_read uses current tags as data size
@ -410,13 +413,13 @@ class ASCIItable():
def info_append(self, def info_append(self,
what): what):
"""Add item or list to existing set of infos""" """Add item or list to existing set of infos"""
if not isinstance(what, (str, unicode)): if isinstance(what, (str, unicode)):
self.info += [self._removeCRLF(what)]
else:
try: try:
for item in what: self.info_append(item) for item in what: self.info_append(item)
except: except TypeError:
self.info += [self._removeCRLF(str(what))] self.info += [self._removeCRLF(str(what))]
else:
self.info += [self._removeCRLF(what)]
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def info_clear(self): def info_clear(self):
@ -468,10 +471,8 @@ class ASCIItable():
"""Read whole data of all (given) labels as numpy array""" """Read whole data of all (given) labels as numpy array"""
from collections import Iterable from collections import Iterable
try: try: self.data_rewind() # try to wind back to start of data
self.data_rewind() # try to wind back to start of data except: pass # assume/hope we are at data start already...
except:
pass # assume/hope we are at data start already...
if labels is None or labels == []: if labels is None or labels == []:
use = None # use all columns (and keep labels intact) use = None # use all columns (and keep labels intact)
@ -530,13 +531,13 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def data_append(self, def data_append(self,
what): what):
if not isinstance(what, (str, unicode)): if isinstance(what, (str, unicode)):
self.data += [what]
else:
try: try:
for item in what: self.data_append(item) for item in what: self.data_append(item)
except: except TypeError:
self.data += [str(what)] self.data += [str(what)]
else:
self.data += [what]
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def data_set(self, def data_set(self,
@ -581,7 +582,7 @@ class ASCIItable():
if len(items) > 2: if len(items) > 2:
if items[1].lower() == 'of': if items[1].lower() == 'of':
items = np.ones(datatype(items[0]))*datatype(items[2]) items = np.ones(datatype(items[0]))*datatype(items[2])
elif items[1].lower() == 'to': elif items[1].lower() == 'to':
items = np.linspace(datatype(items[0]),datatype(items[2]), items = np.linspace(datatype(items[0]),datatype(items[2]),
abs(datatype(items[2])-datatype(items[0]))+1,dtype=int) abs(datatype(items[2])-datatype(items[0]))+1,dtype=int)
else: items = list(map(datatype,items)) else: items = list(map(datatype,items))

2
processing/pre/geom_fromOsteonGeometry.py
View File

@ -65,7 +65,7 @@ if np.any(np.array(options.size) <= 0.0):
if filename == []: filename = [None] if filename == []: filename = [None]
table = damask.ASCIItable(outname = filename[0], table = damask.ASCIItable(outname = filename[0],
buffered = False) buffered = False, labeled=False)
damask.util.report(scriptName,filename[0]) damask.util.report(scriptName,filename[0])

6
src/CMakeLists.txt
View File

@ -17,8 +17,12 @@ list(APPEND OBJECTFILES $<TARGET_OBJECTS:SYSTEM_ROUTINES>)
add_library(PREC OBJECT "prec.f90") add_library(PREC OBJECT "prec.f90")
list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>)
add_library(QUIT OBJECT "quit.f90")
add_dependencies(QUIT PREC)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:QUIT>)
add_library(DAMASK_INTERFACE OBJECT "DAMASK_interface.f90") add_library(DAMASK_INTERFACE OBJECT "DAMASK_interface.f90")
add_dependencies(DAMASK_INTERFACE PREC SYSTEM_ROUTINES) add_dependencies(DAMASK_INTERFACE QUIT SYSTEM_ROUTINES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>)
add_library(IO OBJECT "IO.f90") add_library(IO OBJECT "IO.f90")

407
src/DAMASK_FEM.f90
View File

@ -7,8 +7,13 @@
!> results !> results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
program DAMASK_FEM program DAMASK_FEM
use, intrinsic :: & #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
#include <petsc/finclude/petscsys.h>
use PetscDM
use prec, only: & use prec, only: &
pInt, & pInt, &
pReal, & pReal, &
@ -18,9 +23,7 @@ program DAMASK_FEM
loadCaseFile, & loadCaseFile, &
getSolverJobName getSolverJobName
use IO, only: & use IO, only: &
IO_read, &
IO_isBlank, & IO_isBlank, &
IO_open_file, &
IO_stringPos, & IO_stringPos, &
IO_stringValue, & IO_stringValue, &
IO_floatValue, & IO_floatValue, &
@ -29,12 +32,7 @@ program DAMASK_FEM
IO_lc, & IO_lc, &
IO_intOut, & IO_intOut, &
IO_warning, & IO_warning, &
IO_timeStamp, & IO_timeStamp
IO_EOF
use debug, only: &
debug_level, &
debug_spectral, &
debug_levelBasic
use math ! need to include the whole module for FFTW use math ! need to include the whole module for FFTW
use CPFEM2, only: & use CPFEM2, only: &
CPFEM_initAll CPFEM_initAll
@ -42,9 +40,9 @@ program DAMASK_FEM
restartWrite, & restartWrite, &
restartInc restartInc
use numerics, only: & use numerics, only: &
worldrank, &
maxCutBack, & maxCutBack, &
stagItMax, & stagItMax
worldrank
use mesh, only: & use mesh, only: &
mesh_Nboundaries, & mesh_Nboundaries, &
mesh_boundaries, & mesh_boundaries, &
@ -57,37 +55,17 @@ program DAMASK_FEM
maxFields, & maxFields, &
nActiveFields, & nActiveFields, &
FIELD_MECH_ID, & FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID, &
COMPONENT_MECH_X_ID, & COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, & COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, & FIELD_MECH_label
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID, &
FIELD_MECH_label, &
FIELD_THERMAL_label, &
FIELD_DAMAGE_label, &
FIELD_SOLUTE_label, &
FIELD_MGTWIN_label
use FEM_mech use FEM_mech
implicit none implicit none
#include <petsc/finclude/petsc.h>
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file ! variables related to information from load case and geom file
integer(pInt), parameter :: FILEUNIT = 234_pInt !< file unit, DAMASK IO does not support newunit feature
integer(pInt), allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing integer(pInt), allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing
integer(pInt) :: & integer(pInt) :: &
N_def = 0_pInt !< # of rate of deformation specifiers found in load case file N_def = 0_pInt !< # of rate of deformation specifiers found in load case file
character(len=65536) :: & character(len=65536) :: &
@ -95,7 +73,6 @@ program DAMASK_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! loop variables, convergence etc.
integer(pInt), parameter :: & integer(pInt), parameter :: &
subStepFactor = 2_pInt !< for each substep, divide the last time increment by 2.0 subStepFactor = 2_pInt !< for each substep, divide the last time increment by 2.0
real(pReal) :: & real(pReal) :: &
@ -105,7 +82,8 @@ program DAMASK_FEM
timeIncOld = 0.0_pReal, & !< previous time interval timeIncOld = 0.0_pReal, & !< previous time interval
remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case
logical :: & logical :: &
guess !< guess along former trajectory guess, & !< guess along former trajectory
stagIterate
integer(pInt) :: & integer(pInt) :: &
i, & i, &
errorID, & errorID, &
@ -115,18 +93,18 @@ program DAMASK_FEM
currentLoadcase = 0_pInt, & !< current load case currentLoadcase = 0_pInt, & !< current load case
currentFace = 0_pInt, & currentFace = 0_pInt, &
inc, & !< current increment in current load case inc, & !< current increment in current load case
totalIncsCounter = 0_pInt, & !< total No. of increments totalIncsCounter = 0_pInt, & !< total # of increments
convergedCounter = 0_pInt, & !< No. of converged increments convergedCounter = 0_pInt, & !< # of converged increments
notConvergedCounter = 0_pInt, & !< No. of non-converged increments notConvergedCounter = 0_pInt, & !< # of non-converged increments
fileUnit = 0_pInt, & !< file unit for reading load case and writing results
myStat, &
statUnit = 0_pInt, & !< file unit for statistics output statUnit = 0_pInt, & !< file unit for statistics output
lastRestartWritten = 0_pInt !< total increment No. at which last restart information was written lastRestartWritten = 0_pInt, & !< total increment No. at which last restart information was written
integer(pInt) :: &
stagIter, & stagIter, &
component component
logical :: &
stagIterate
character(len=6) :: loadcase_string character(len=6) :: loadcase_string
character(len=1024) :: incInfo !< string parsed to solution with information about current load case character(len=1024) :: &
incInfo
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet PetscInt :: faceSet, currentFaceSet
@ -134,17 +112,13 @@ program DAMASK_FEM
PetscErrorCode :: ierr PetscErrorCode :: ierr
external :: & external :: &
MPI_abort, &
DMGetDimension, &
DMGetLabelSize, &
DMGetLabelIdIS, &
ISDestroy, &
quit quit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll call CPFEM_initAll
write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
! reading basic information from load case file and allocate data structure containing load cases ! reading basic information from load case file and allocate data structure containing load cases
@ -154,12 +128,13 @@ program DAMASK_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases ! reading basic information from load case file and allocate data structure containing load cases
call IO_open_file(FILEUNIT,trim(loadCaseFile)) open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read')
rewind(FILEUNIT) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=trim(loadCaseFile))
do do
line = IO_read(FILEUNIT) read(fileUnit, '(A)', iostat=myStat) line
if (trim(line) == IO_EOF) exit if ( myStat /= 0_pInt) exit
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1) ! reading compulsory parameters for loadcase do i = 1_pInt, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i))) select case (IO_lc(IO_stringValue(line,chunkPos,i)))
@ -183,6 +158,16 @@ program DAMASK_FEM
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements
allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents)) allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents))
do component = 1, loadCases(i)%fieldBC(field)%nComponents
select case (component)
case (1)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_X_ID
case (2)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
case (3)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select
enddo
end select end select
do component = 1, loadCases(i)%fieldBC(field)%nComponents do component = 1, loadCases(i)%fieldBC(field)%nComponents
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal) allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
@ -193,11 +178,12 @@ program DAMASK_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading the load case and assign values to the allocated data structure ! reading the load case and assign values to the allocated data structure
rewind(FILEUNIT) rewind(fileUnit)
do do
line = IO_read(FILEUNIT) read(fileUnit, '(A)', iostat=myStat) line
if (trim(line) == IO_EOF) exit if ( myStat /= 0_pInt) exit
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1) do i = 1_pInt, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i))) select case (IO_lc(IO_stringValue(line,chunkPos,i)))
@ -268,134 +254,15 @@ program DAMASK_FEM
enddo enddo
endif endif
enddo enddo
case('temp','temperature') ! thermal field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_THERMAL_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_THERMAL_T_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('mgtwin') ! mgtwin field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MGTWIN_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MGTWIN_PHI_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('damage')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_DAMAGE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_DAMAGE_PHI_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cv')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CV_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('ch')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CH_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('chpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CHPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvah')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaH_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvahpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaHPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
end select end select
enddo; enddo enddo; enddo
close(FILEUNIT) close(fileUnit)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! consistency checks and output of load case ! consistency checks and output of load case
loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase
errorID = 0_pInt errorID = 0_pInt
if (worldrank == 0) then checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases)
checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases)
write (loadcase_string, '(i6)' ) currentLoadCase write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
@ -405,18 +272,6 @@ program DAMASK_FEM
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label) write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label)
case(FIELD_THERMAL_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_THERMAL_label)
case(FIELD_DAMAGE_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_DAMAGE_label)
case(FIELD_MGTWIN_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MGTWIN_label)
case(FIELD_SOLUTE_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_SOLUTE_label)
end select end select
do faceSet = 1_pInt, mesh_Nboundaries do faceSet = 1_pInt, mesh_Nboundaries
do component = 1_pInt, loadCases(currentLoadCase)%fieldBC(field)%nComponents do component = 1_pInt, loadCases(currentLoadCase)%fieldBC(field)%nComponents
@ -437,11 +292,10 @@ program DAMASK_FEM
write(6,'(2x,a,i5,/)') 'restart frequency: ', & write(6,'(2x,a,i5,/)') 'restart frequency: ', &
loadCases(currentLoadCase)%restartfrequency loadCases(currentLoadCase)%restartfrequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases enddo checkLoadcases
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing initialization depending on selected solver ! doing initialization depending on active solvers
call Utilities_init() call Utilities_init()
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(1)%fieldBC(field)%ID) select case (loadCases(1)%fieldBC(field)%ID)
@ -450,56 +304,49 @@ program DAMASK_FEM
end select end select
enddo enddo
!--------------------------------------------------------------------------------------------------
! loopping over loadcases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
time0 = time ! currentLoadCase start time
if (loadCases(currentLoadCase)%followFormerTrajectory) then
guess = .true.
else
guess = .false. ! change of load case, homogeneous guess for the first inc
endif
!-------------------------------------------------------------------------------------------------- loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
! loop oper incs defined in input file for current currentLoadCase time0 = time ! load case start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt totalIncsCounter = totalIncsCounter + 1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forwarding time ! forwarding time
timeIncOld = timeinc timeIncOld = timeinc ! last timeinc that brought former inc to an end
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/loadCases(currentLoadCase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else else
if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale if (currentLoadCase == 1_pInt) then ! 1st load case of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale if (inc == 1_pInt) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st currentLoadCase of logarithmic scale else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal)) timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal))
endif endif
else ! not-1st currentLoadCase of logarithmic scale else ! not-1st load case of logarithmic scale
timeinc = time0 * & timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/& ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))& real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/& -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1_pInt ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))) real(loadCases(currentLoadCase)%incs ,pReal)))
endif endif
endif endif
timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
forwarding: if(totalIncsCounter >= restartInc) then skipping: if (totalIncsCounter <= restartInc) then ! not yet at restart inc?
stepFraction = 0_pInt time = time + timeinc ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping
stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel
!-------------------------------------------------------------------------------------------------- subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
! loop over sub incs remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt) time = time + timeinc ! forward target time
time = time + timeinc ! forward time stepFraction = stepFraction + 1_pInt ! count step
stepFraction = stepFraction + 1_pInt
remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! report begin of new step
if (worldrank == 0) then
write(6,'(/,a)') ' ###########################################################################' write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//& write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
@ -509,12 +356,14 @@ program DAMASK_FEM
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,& '-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases) ' of load case ', currentLoadCase,'/',size(loadCases)
flush(6) write(incInfo,&
write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//& '(a,'//IO_intOut(totalIncsCounter)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') & ',a,'//IO_intOut(sum(loadCases%incs))//&
',a,'//IO_intOut(stepFraction)//&
',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel '-',stepFraction, '/', subStepFactor**cutBackLevel
endif flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forward fields ! forward fields
@ -539,12 +388,14 @@ program DAMASK_FEM
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select end select
if(.not. solres(field)%converged) exit ! no solution found if(.not. solres(field)%converged) exit ! no solution found
enddo enddo
stagIter = stagIter + 1_pInt stagIter = stagIter + 1_pInt
stagIterate = stagIter < stagItMax .and. & stagIterate = stagIter < stagItMax &
all(solres(:)%converged) .and. & .and. all(solres(:)%converged) &
.not. all(solres(:)%stagConverged) .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
enddo enddo
! check solution ! check solution
@ -570,85 +421,49 @@ program DAMASK_FEM
if (worldrank == 0) write(statUnit,*) totalIncsCounter, time, cutBackLevel, & if (worldrank == 0) write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
endif endif
enddo subIncLooping enddo subStepLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
if(all(solres(:)%converged)) then ! report converged inc
if (all(solres(:)%converged)) then
convergedCounter = convergedCounter + 1_pInt convergedCounter = convergedCounter + 1_pInt
if (worldrank == 0) then write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ' increment ', totalIncsCounter, ' converged'
' increment ', totalIncsCounter, ' converged'
endif
else else
if (worldrank == 0) then
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif
notConvergedCounter = notConvergedCounter + 1_pInt notConvergedCounter = notConvergedCounter + 1_pInt
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6) endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
if (worldrank == 0) then write(6,'(1/,a)') ' ToDo: ... writing results to file ......................................'
write(6,'(1/,a)') ' ... writing results to file ......................................'
endif
endif endif
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ...
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ToDo first call to CPFEM_general will write? .and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information
restartWrite = .true. restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif endif
else forwarding
time = time + timeinc endif skipping
guess = .true.
endif forwarding
enddo incLooping enddo incLooping
enddo loadCaseLooping enddo loadCaseLooping
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report summary of whole calculation ! report summary of whole calculation
if (worldrank == 0) then
write(6,'(/,a)') ' ###########################################################################' write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & write(6,'(1x,'//IO_intOut(convergedCounter)//',a,'//IO_intOut(notConvergedCounter + convergedCounter)//',a,f5.1,a)') &
notConvergedCounter + convergedCounter, ' (', & convergedCounter, ' out of ', &
real(convergedCounter, pReal)/& notConvergedCounter + convergedCounter, ' (', &
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & real(convergedCounter, pReal)/&
' %) increments converged!' real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!'
endif flush(6)
if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged call MPI_file_close(fileUnit,ierr)
call quit(0_pInt) ! no complains ;) close(statUnit)
end program DAMASK_FEM if (notConvergedCounter > 0_pInt) call quit(2_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;)
end program DAMASK_FEM
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief quit subroutine to mimic behavior of FEM solvers
!> @details exits the Spectral solver and reports time and duration. Exit code 0 signals
!> everything went fine. Exit code 1 signals an error, message according to IO_error. Exit code
!> 2 signals request for regridding, increment of last saved restart information is written to
!> stderr. Exit code 3 signals no severe problems, but some increments did not converge
!--------------------------------------------------------------------------------------------------
subroutine quit(stop_id)
use prec, only: &
pInt
implicit none
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id < 0_pInt) then ! trigger regridding
write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
stop 2
endif
if (stop_id == 3_pInt) stop 3 ! not all incs converged
stop 1 ! error (message from IO_error)
end subroutine quit

50
src/DAMASK_interface.f90
View File

@ -43,24 +43,23 @@ subroutine DAMASK_interface_init()
use, intrinsic :: & use, intrinsic :: &
iso_fortran_env iso_fortran_env
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
#if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR!=9 #if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR!=10
=================================================================================================== ===================================================================================================
3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x
=================================================================================================== ===================================================================================================
======= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =========================================== ======= THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ==========================================
========== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ======================================== ========== THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x =======================================
============= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ===================================== ============= THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ====================================
================ THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ================================== ================ THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x =================================
=================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =============================== =================== THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ==============================
====================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ============================ ====================== THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ===========================
========================= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ========================= ========================= THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ========================
============================ THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ====================== ============================ THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x =====================
=============================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =================== =============================== THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ==================
================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ================ ================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ===============
===================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ============= ===================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.10.x ============
======================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ==========
=================================================================================================== ===================================================================================================
3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x 3.10.x
=================================================================================================== ===================================================================================================
#endif #endif
use PETScSys use PETScSys
@ -88,9 +87,7 @@ subroutine DAMASK_interface_init()
dateAndTime ! type default integer dateAndTime ! type default integer
PetscErrorCode :: ierr PetscErrorCode :: ierr
external :: & external :: &
quit,& quit
PETScErrorF, & ! is called in the CHKERRQ macro
PETScInitialize
open(6, encoding='UTF-8') ! for special characters in output open(6, encoding='UTF-8') ! for special characters in output
@ -230,8 +227,8 @@ subroutine setWorkingDirectory(workingDirectoryArg)
implicit none implicit none
character(len=*), intent(in) :: workingDirectoryArg !< working directory argument character(len=*), intent(in) :: workingDirectoryArg !< working directory argument
character(len=1024) :: workingDirectory !< working directory argument character(len=1024) :: workingDirectory !< working directory argument
external :: quit
logical :: error logical :: error
external :: quit
absolutePath: if (workingDirectoryArg(1:1) == '/') then absolutePath: if (workingDirectoryArg(1:1) == '/') then
workingDirectory = workingDirectoryArg workingDirectory = workingDirectoryArg
@ -284,11 +281,18 @@ character(len=1024) function getGeometryFile(geometryParameter)
implicit none implicit none
character(len=1024), intent(in) :: geometryParameter character(len=1024), intent(in) :: geometryParameter
logical :: file_exists
external :: quit
getGeometryFile = trim(geometryParameter) getGeometryFile = trim(geometryParameter)
if (scan(getGeometryFile,'/') /= 1) getGeometryFile = trim(getCWD())//'/'//trim(getGeometryFile) if (scan(getGeometryFile,'/') /= 1) getGeometryFile = trim(getCWD())//'/'//trim(getGeometryFile)
getGeometryFile = makeRelativePath(trim(getCWD()), getGeometryFile) getGeometryFile = makeRelativePath(trim(getCWD()), getGeometryFile)
inquire(file=trim(getGeometryFile), exist=file_exists)
if (.not. file_exists) then
write(6,'(a)') ' Geometry file does not exists ('//trim(getGeometryFile)//')'
call quit(1_pInt)
endif
end function getGeometryFile end function getGeometryFile
@ -302,11 +306,19 @@ character(len=1024) function getLoadCaseFile(loadCaseParameter)
implicit none implicit none
character(len=1024), intent(in) :: loadCaseParameter character(len=1024), intent(in) :: loadCaseParameter
logical :: file_exists
external :: quit
getLoadCaseFile = trim(loadCaseParameter) getLoadCaseFile = trim(loadCaseParameter)
if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = trim(getCWD())//'/'//trim(getLoadCaseFile) if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = trim(getCWD())//'/'//trim(getLoadCaseFile)
getLoadCaseFile = makeRelativePath(trim(getCWD()), getLoadCaseFile) getLoadCaseFile = makeRelativePath(trim(getCWD()), getLoadCaseFile)
inquire(file=trim(getLoadCaseFile), exist=file_exists)
if (.not. file_exists) then
write(6,'(a)') ' Geometry file does not exists ('//trim(getLoadCaseFile)//')'
call quit(1_pInt)
endif
end function getLoadCaseFile end function getLoadCaseFile

75
src/DAMASK_spectral.f90
View File

@ -88,7 +88,6 @@ program DAMASK_spectral
real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0) real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: & integer(pInt) :: &
N_t = 0_pInt, & !< # of time indicators found in load case file N_t = 0_pInt, & !< # of time indicators found in load case file
N_n = 0_pInt, & !< # of increment specifiers found in load case file N_n = 0_pInt, & !< # of increment specifiers found in load case file
@ -130,8 +129,7 @@ program DAMASK_spectral
stagIter stagIter
character(len=6) :: loadcase_string character(len=6) :: loadcase_string
character(len=1024) :: & character(len=1024) :: &
incInfo, & !< string parsed to solution with information about current load case incInfo
workingDir
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tLoadCase) :: newLoadCase type(tLoadCase) :: newLoadCase
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
@ -140,7 +138,7 @@ program DAMASK_spectral
integer(pInt), parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742 integer(pInt), parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742
integer(pInt), parameter :: maxRealOut = maxByteOut/pReal integer(pInt), parameter :: maxRealOut = maxByteOut/pReal
integer(pLongInt), dimension(2) :: outputIndex integer(pLongInt), dimension(2) :: outputIndex
integer :: ierr PetscErrorCode :: ierr
procedure(basic_init), pointer :: & procedure(basic_init), pointer :: &
mech_init mech_init
procedure(basic_forward), pointer :: & procedure(basic_forward), pointer :: &
@ -151,7 +149,6 @@ program DAMASK_spectral
external :: & external :: &
quit quit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll call CPFEM_initAll
@ -378,7 +375,7 @@ program DAMASK_spectral
open(newunit=fileUnit,file=trim(getSolverJobName())//& open(newunit=fileUnit,file=trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED',status='REPLACE') '.spectralOut',form='UNFORMATTED',status='REPLACE')
write(fileUnit) 'load:', trim(loadCaseFile) ! ... and write header write(fileUnit) 'load:', trim(loadCaseFile) ! ... and write header
write(fileUnit) 'workingdir:', trim(workingDir) write(fileUnit) 'workingdir:', 'n/a'
write(fileUnit) 'geometry:', trim(geometryFile) write(fileUnit) 'geometry:', trim(geometryFile)
write(fileUnit) 'grid:', grid write(fileUnit) 'grid:', grid
write(fileUnit) 'size:', geomSize write(fileUnit) 'size:', geomSize
@ -434,14 +431,12 @@ program DAMASK_spectral
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif writeUndeformed endif writeUndeformed
!--------------------------------------------------------------------------------------------------
! looping over load cases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
time0 = time ! load case start time time0 = time ! load case start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc
!--------------------------------------------------------------------------------------------------
! loop over incs defined in input file for current load case
incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt totalIncsCounter = totalIncsCounter + 1_pInt
@ -473,8 +468,6 @@ program DAMASK_spectral
else skipping else skipping
stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel
!--------------------------------------------------------------------------------------------------
! loop over sub step
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time time = time + timeinc ! forward target time
@ -629,63 +622,7 @@ program DAMASK_spectral
call MPI_file_close(fileUnit,ierr) call MPI_file_close(fileUnit,ierr)
close(statUnit) close(statUnit)
if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged if (notConvergedCounter > 0_pInt) call quit(2_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;) call quit(0_pInt) ! no complains ;)
end program DAMASK_spectral end program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief quit subroutine to mimic behavior of FEM solvers
!> @details exits the Spectral solver and reports time and duration. Exit code 0 signals
!> everything went fine. Exit code 1 signals an error, message according to IO_error. Exit code
!> 2 signals no converged solution and increment of last saved restart information is written to
!> stderr. Exit code 3 signals no severe problems, but some increments did not converge
!--------------------------------------------------------------------------------------------------
subroutine quit(stop_id)
#include <petsc/finclude/petscsys.h>
#ifdef _OPENMP
use MPI, only: &
MPI_finalize
#endif
use prec, only: &
pInt
implicit none
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
integer(pInt) :: error = 0_pInt
PetscErrorCode :: ierr = 0
logical :: ErrorInQuit
external :: &
PETScFinalize
call PETScFinalize(ierr)
if (ierr /= 0) write(6,'(a)') ' Error in PETScFinalize'
#ifdef _OPENMP
call MPI_finalize(error)
if (error /= 0) write(6,'(a)') ' Error in MPI_finalize'
#endif
ErrorInQuit = (ierr /= 0 .or. error /= 0_pInt)
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
if (stop_id == 0_pInt .and. .not. ErrorInQuit) stop 0 ! normal termination
if (stop_id < 0_pInt .and. .not. ErrorInQuit) then ! terminally ill, restart might help
write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
stop 2
endif
if (stop_id == 3_pInt .and. .not. ErrorInQuit) stop 3 ! not all incs converged
stop 1 ! error (message from IO_error)
end subroutine quit

89
src/FEM_mech.f90
View File

@ -5,12 +5,14 @@
!> @brief FEM PETSc solver !> @brief FEM PETSc solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module FEM_mech module FEM_mech
#include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscdm.h>
#include <petsc/finclude/petsc.h> #include <petsc/finclude/petsc.h>
use PETScdmda use PETScsnes
use PETScsnes use PETScDM
use PETScDM use PETScDMplex
use PETScDMplex use PETScDT
use prec, only: & use prec, only: &
pInt, & pInt, &
pReal pReal
@ -20,9 +22,6 @@ use PETScDMplex
tSolutionState, & tSolutionState, &
tFieldBC, & tFieldBC, &
tComponentBC tComponentBC
use numerics, only: &
worldrank, &
worldsize
use mesh, only: & use mesh, only: &
mesh_Nboundaries, & mesh_Nboundaries, &
mesh_boundaries mesh_boundaries
@ -45,7 +44,7 @@ use PETScDMplex
SNES, private :: mech_snes SNES, private :: mech_snes
Vec, private :: solution, solution_rate, solution_local Vec, private :: solution, solution_rate, solution_local
PetscInt, private :: dimPlex, cellDof, nQuadrature, nBasis PetscInt, private :: dimPlex, cellDof, nQuadrature, nBasis
PetscReal, allocatable, target,dimension(:), private :: qPoints, qWeights PetscReal, allocatable, target, private :: qPoints(:), qWeights(:)
MatNullSpace, private :: matnull MatNullSpace, private :: matnull
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -61,27 +60,6 @@ use PETScDMplex
FEM_mech_solution ,& FEM_mech_solution ,&
FEM_mech_forward, & FEM_mech_forward, &
FEM_mech_destroy FEM_mech_destroy
external :: &
MatZeroRowsColumnsLocalIS, &
PetscQuadratureCreate, &
PetscFECreateDefault, &
PetscFESetQuadrature, &
PetscFEGetDimension, &
PetscFEDestroy, &
PetscFEGetDualSpace, &
PetscQuadratureDestroy, &
PetscDSSetDiscretization, &
PetscDSGetTotalDimension, &
PetscDSGetDiscretization, &
PetscDualSpaceGetFunctional, &
DMGetLabelSize, &
DMSNESSetFunctionLocal, &
DMSNESSetJacobianLocal, &
SNESSetOptionsPrefix, &
SNESSetConvergenceTest, &
PetscObjectSetName
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -97,7 +75,6 @@ subroutine FEM_mech_init(fieldBC)
use mesh, only: & use mesh, only: &
geomMesh geomMesh
use numerics, only: & use numerics, only: &
worldrank, &
itmax, & itmax, &
integrationOrder integrationOrder
use FEM_Zoo, only: & use FEM_Zoo, only: &
@ -115,14 +92,14 @@ subroutine FEM_mech_init(fieldBC)
DMLabel :: BCLabel DMLabel :: BCLabel
PetscInt, allocatable, target :: numComp(:), numDoF(:), bcField(:) PetscInt, allocatable, target :: numComp(:), numDoF(:), bcField(:)
PetscInt, pointer :: pNumComp(:), pNumDof(:), pBcField(:), pBcPoint(:) PetscInt, pointer :: pNumComp(:), pNumDof(:), pBcField(:), pBcPoint(:)
PetscInt :: numBC, bcSize PetscInt :: numBC, bcSize, nc
IS :: bcPoint IS :: bcPoint
IS, allocatable, target :: bcComps(:), bcPoints(:) IS, allocatable, target :: bcComps(:), bcPoints(:)
IS, pointer :: pBcComps(:), pBcPoints(:) IS, pointer :: pBcComps(:), pBcPoints(:)
PetscSection :: section PetscSection :: section
PetscInt :: field, faceSet, topologDim, nNodalPoints PetscInt :: field, faceSet, topologDim, nNodalPoints
PetscReal, pointer :: qPointsP(:), qWeightsP(:), & PetscReal, dimension(:) , pointer :: qPointsP, qWeightsP, &
nodalPointsP(:), nodalWeightsP(:) nodalPointsP, nodalWeightsP
PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:) PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:)
PetscScalar, pointer :: px_scal(:) PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:) PetscScalar, allocatable, target :: x_scal(:)
@ -132,7 +109,7 @@ subroutine FEM_mech_init(fieldBC)
PetscInt :: cellStart, cellEnd, cell, basis PetscInt :: cellStart, cellEnd, cell, basis
character(len=7) :: prefix = 'mechFE_' character(len=7) :: prefix = 'mechFE_'
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>' write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
@ -144,20 +121,21 @@ subroutine FEM_mech_init(fieldBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech discretization ! Setup FEM mech discretization
allocate(qPoints(dimPlex*FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
allocate(qWeights(FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
qPoints = FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p qPoints = FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p
qWeights = FEM_Zoo_QuadratureWeights(dimPlex,integrationOrder)%p qWeights = FEM_Zoo_QuadratureWeights(dimPlex,integrationOrder)%p
nQuadrature = FEM_Zoo_nQuadrature(dimPlex,integrationOrder) nQuadrature = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
qPointsP => qPoints qPointsP => qPoints
qWeightsP => qWeights qWeightsP => qWeights
call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr) call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr)
call PetscQuadratureSetData(mechQuad,dimPlex,nQuadrature,qPointsP,qWeightsP,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscFECreateDefault(mech_mesh,dimPlex,dimPlex,PETSC_TRUE,prefix, & nc = dimPlex
call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr)
CHKERRQ(ierr)
call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, &
integrationOrder,mechFE,ierr); CHKERRQ(ierr) integrationOrder,mechFE,ierr); CHKERRQ(ierr)
call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
nBasis = nBasis/nc
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr) call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr)
call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr) call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr)
@ -168,7 +146,7 @@ subroutine FEM_mech_init(fieldBC)
! Setup FEM mech boundary conditions ! Setup FEM mech boundary conditions
call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr) call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr)
call DMGetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
allocate(numComp(1), source=dimPlex); pNumComp => numComp allocate(numComp(1), source=dimPlex); pNumComp => numComp
allocate(numDof(dimPlex+1), source = 0); pNumDof => numDof allocate(numDof(dimPlex+1), source = 0); pNumDof => numDof
do topologDim = 0, dimPlex do topologDim = 0, dimPlex
@ -210,7 +188,7 @@ subroutine FEM_mech_init(fieldBC)
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS, & numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS, &
section,ierr) section,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
do faceSet = 1, numBC do faceSet = 1, numBC
call ISDestroy(bcPoints(faceSet),ierr); CHKERRQ(ierr) call ISDestroy(bcPoints(faceSet),ierr); CHKERRQ(ierr)
enddo enddo
@ -250,7 +228,7 @@ subroutine FEM_mech_init(fieldBC)
pV0 => v0 pV0 => v0
pCellJ => cellJ pCellJ => cellJ
pInvcellJ => invcellJ pInvcellJ => invcellJ
call DMGetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSGetDiscretization(mechDS,0,mechFE,ierr) call PetscDSGetDiscretization(mechDS,0,mechFE,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -262,12 +240,12 @@ subroutine FEM_mech_init(fieldBC)
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
do basis = 0, nBasis-1 do basis = 0, nBasis*dimPlex-1, dimPlex
call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr) call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscQuadratureGetData(functional,dimPlex,nNodalPoints,nodalPointsP,nodalWeightsP,ierr) call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
x_scal(basis*dimPlex+1:(basis+1)*dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0) x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0)
enddo enddo
px_scal => x_scal px_scal => x_scal
call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr) call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr)
@ -371,7 +349,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
pV0 => v0 pV0 => v0
pCellJ => cellJ pCellJ => cellJ
pInvcellJ => invcellJ pInvcellJ => invcellJ
call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr) call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr) call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -405,8 +383,8 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx )*dimPlex+1: & matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: &
(qPt*nBasis*dimPlex+cidx+1)*dimPlex )) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
@ -446,8 +424,8 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx )*dimPlex+1: & matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: &
(qPt*nBasis*dimPlex+cidx+1)*dimPlex )) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
f_scal = f_scal + & f_scal = f_scal + &
@ -520,8 +498,8 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) call MatZeroEntries(Jac,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr) call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr) call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetDefaultGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr) call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr) call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
@ -555,8 +533,8 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx )*dimPlex+1: & matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: &
(qPt*nBasis*dimPlex+cidx+1)*dimPlex )) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), &
@ -605,7 +583,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! apply boundary conditions ! apply boundary conditions
!call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr) MD: linker error call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr)
call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr) call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)
@ -646,7 +624,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
ForwardData = .True. ForwardData = .True.
materialpoint_F0 = materialpoint_F materialpoint_F0 = materialpoint_F
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local
call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0,ierr); CHKERRQ(ierr) call VecSet(x_local,0.0,ierr); CHKERRQ(ierr)
call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector
@ -684,6 +662,7 @@ end subroutine FEM_mech_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: & use numerics, only: &
worldrank, &
err_struct_tolAbs, & err_struct_tolAbs, &
err_struct_tolRel err_struct_tolRel
use IO, only: & use IO, only: &

240
src/FEM_utilities.f90
View File

@ -3,9 +3,12 @@
!> @brief Utilities used by the FEM solver !> @brief Utilities used by the FEM solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module FEM_utilities module FEM_utilities
#include <petsc/finclude/petsc.h> #include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscdmda.h>
#include <petsc/finclude/petscis.h>
use prec, only: pReal, pInt use prec, only: pReal, pInt
use PETScdmplex
use PETScdmda use PETScdmda
use PETScis use PETScis
@ -21,25 +24,15 @@ use PETScis
! grid related information information ! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems real(pReal), public :: wgt !< weighting factor 1/Nelems
real(pReal), public :: wgtDof !< weighting factor 1/Nelems real(pReal), public :: wgtDof !< weighting factor 1/Nelems
real(pReal), public :: C_volAvg(3,3,3,3)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output data ! output data
PetscViewer, public :: resUnit
Vec, public :: coordinatesVec Vec, public :: coordinatesVec
Vec, allocatable, public :: homogenizationResultsVec(:), &
crystalliteResultsVec(:,:), &
phaseResultsVec(:,:)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! field labels information ! field labels information
character(len=*), parameter, public :: & character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical', & FIELD_MECH_label = 'mechanical'
FIELD_THERMAL_label = 'thermal', &
FIELD_DAMAGE_label = 'damage', &
FIELD_SOLUTE_label = 'solute', &
FIELD_MGTWIN_label = 'mgtwin'
enum, bind(c) enum, bind(c)
enumerator :: FIELD_UNDEFINED_ID, & enumerator :: FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, & FIELD_MECH_ID, &
@ -123,34 +116,10 @@ use PETScis
utilities_indexActiveSet, & utilities_indexActiveSet, &
utilities_destroy, & utilities_destroy, &
FIELD_MECH_ID, & FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID, &
COMPONENT_MECH_X_ID, & COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, & COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, & COMPONENT_THERMAL_T_ID
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID
external :: &
PetscOptionsInsertString, &
PetscObjectSetName, &
DMPlexGetHeightStratum, &
DMGetLabelIdIS, &
DMPlexGetChart, &
DMPlexLabelComplete, &
PetscViewerHDF5Open, &
PetscViewerHDF5PushGroup, &
PetscViewerHDF5PopGroup, &
PetscViewerDestroy
contains contains
@ -167,6 +136,7 @@ subroutine utilities_init()
IO_timeStamp, & IO_timeStamp, &
IO_open_file IO_open_file
use numerics, only: & use numerics, only: &
structOrder, &
integrationOrder, & integrationOrder, &
worldsize, & worldsize, &
worldrank, & worldrank, &
@ -190,14 +160,11 @@ subroutine utilities_init()
implicit none implicit none
character(len=1024) :: petsc_optionsPhysics, grainStr character(len=1024) :: petsc_optionsPhysics
integer(pInt) :: dimPlex integer(pInt) :: dimPlex
integer(pInt) :: headerID = 205_pInt integer(pInt) :: headerID = 205_pInt
PetscInt, dimension(:), pointer :: points PetscInt, allocatable :: nEntities(:), nOutputCells(:), nOutputNodes(:)
PetscInt, allocatable :: nEntities(:), nOutputCells(:), nOutputNodes(:), mappingCells(:) PetscInt :: dim
PetscInt :: cellStart, cellEnd, cell, ip, dim, ctr, qPt
PetscInt, allocatable :: connectivity(:,:)
Vec :: connectivityVec
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
@ -221,15 +188,12 @@ subroutine utilities_init()
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
!write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order ' , structOrder write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_degree ' , structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal) wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
call PetscViewerHDF5Open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.h5', &
FILE_MODE_WRITE, resUnit, ierr); CHKERRQ(ierr)
call PetscViewerHDF5PushGroup(resUnit, '/', ierr); CHKERRQ(ierr)
call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRQ(ierr) call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRQ(ierr)
allocate(nEntities(dimPlex+1), source=0) allocate(nEntities(dimPlex+1), source=0)
allocate(nOutputNodes(worldsize), source = 0) allocate(nOutputNodes(worldsize), source = 0)
@ -257,178 +221,6 @@ subroutine utilities_init()
write(headerID, '(a,i0)') 'number of cells : ', sum(nOutputCells) write(headerID, '(a,i0)') 'number of cells : ', sum(nOutputCells)
endif endif
allocate(connectivity(2**dimPlex,nOutputCells(worldrank+1)))
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr)
ctr = 0
select case (integrationOrder)
case(1_pInt)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexGetTransitiveClosure(geomMesh,cell,PETSC_TRUE,points,ierr)
CHKERRQ(ierr)
if (dimPlex == 2) then
connectivity(:,ctr+1) = [points( 9), points(11), points(13), points(13)] - nEntities(dimPlex+1)
ctr = ctr + 1
else
connectivity(:,ctr+1) = [points(23), points(25), points(27), points(27), &
points(29), points(29), points(29), points(29)] - nEntities(dimPlex+1)
ctr = ctr + 1
endif
enddo
case(2_pInt)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexGetTransitiveClosure(geomMesh,cell,PETSC_TRUE,points,ierr)
CHKERRQ(ierr)
if (dimPlex == 2) then
connectivity(:,ctr+1) = [points(9 ), points(3), points(1), points(7)]
connectivity(:,ctr+2) = [points(11), points(5), points(1), points(3)]
connectivity(:,ctr+3) = [points(13), points(7), points(1), points(5)]
ctr = ctr + 3
else
connectivity(:,ctr+1) = [points(23), points(11), points(3), points(15), points(17), points(5), points(1), points(7)]
connectivity(:,ctr+2) = [points(25), points(13), points(3), points(11), points(19), points(9), points(1), points(5)]
connectivity(:,ctr+3) = [points(27), points(15), points(3), points(13), points(21), points(7), points(1), points(9)]
connectivity(:,ctr+4) = [points(29), points(17), points(7), points(21), points(19), points(5), points(1), points(9)]
ctr = ctr + 4_pInt
endif
enddo
case default
do cell = cellStart, cellEnd-1; do ip = 0, mesh_maxNips-1
connectivity(:,ctr+1) = cell*mesh_maxNips + ip
ctr = ctr + 1
enddo; enddo
end select
connectivity = connectivity + sum(nOutputNodes(1:worldrank))
call VecCreateMPI(PETSC_COMM_WORLD,dimPlex*nOutputNodes(worldrank+1),dimPlex*sum(nOutputNodes), &
coordinatesVec,ierr);CHKERRQ(ierr)
call PetscObjectSetName(coordinatesVec, 'NodalCoordinates',ierr)
call VecSetFromOptions(coordinatesVec, ierr); CHKERRQ(ierr)
!allocate(mappingCells(worldsize), source = 0)
!do homog = 1, material_Nhomogenization
! mappingCells = 0_pInt; mappingCells(worldrank+1) = homogOutput(homog)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! call PetscObjectSetName(connectivityVec,'mapping_'//trim(homogenization_name(homog)),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec,results,ierr); CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (material_homog(qPt,cell+1) == homog) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo
!do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
! mappingCells = 0_pInt
! mappingCells(worldrank+1) = crystalliteOutput(cryst,grain)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! call PetscObjectSetName(connectivityVec,'mapping_'// &
! trim(crystallite_name(cryst))//'_'// &
! trim(grainStr),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (homogenization_Ngrains (mesh_element(3,cell+1)) >= grain .and. &
! microstructure_crystallite(mesh_element(4,cell+1)) == cryst) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo; enddo
!do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
! mappingCells = 0_pInt
! mappingCells(worldrank+1) = phaseOutput(phase,grain)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! call PetscObjectSetName(connectivityVec,&
! 'mapping_'//trim(phase_name(phase))//'_'// &
! trim(grainStr),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec, results, ierr)
! CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (material_phase(grain,qPt,cell+1) == phase) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr)
! CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr);CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr);CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr);CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo; enddo
!if (worldrank == 0_pInt) then
! do homog = 1, material_Nhomogenization
! call VecGetSize(homogenizationResultsVec(homog),mappingCells(1),ierr)
! CHKERRQ(ierr)
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of homog_'// &
! trim(homogenization_name(homog))//'_'// &
! 'cells : ', mappingCells(1)
! enddo
! do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
! call VecGetSize(crystalliteResultsVec(cryst,grain),mappingCells(1),ierr)
! CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of cryst_'// &
! trim(crystallite_name(cryst))//'_'// &
! trim(grainStr)//'_'// &
! 'cells : ', mappingCells(1)
! enddo; enddo
! do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
! call VecGetSize(phaseResultsVec(phase,grain),mappingCells(1),ierr)
! CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of phase_'// &
! trim(phase_name(phase))//'_'//trim(grainStr)//'_'// &
! 'cells : ', mappingCells(1)
! enddo; enddo
! close(headerID)
!endif
end subroutine utilities_init end subroutine utilities_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -438,8 +230,6 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
use debug, only: & use debug, only: &
debug_reset, & debug_reset, &
debug_info debug_info
use numerics, only: &
worldrank
use math, only: & use math, only: &
math_transpose33, & math_transpose33, &
math_rotate_forward33, & math_rotate_forward33, &
@ -447,10 +237,8 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
use FEsolving, only: & use FEsolving, only: &
restartWrite restartWrite
use homogenization, only: & use homogenization, only: &
materialpoint_F0, &
materialpoint_F, & materialpoint_F, &
materialpoint_P, & materialpoint_P, &
materialpoint_dPdF, &
materialpoint_stressAndItsTangent materialpoint_stressAndItsTangent
use mesh, only: & use mesh, only: &
mesh_NcpElems mesh_NcpElems
@ -503,9 +291,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
cutBack = .false. ! reset cutBack status cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD, ierr)
end subroutine utilities_constitutiveResponse end subroutine utilities_constitutiveResponse

57
src/constitutive.f90
View File

@ -440,7 +440,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
math_Mandel33to6, & math_Mandel33to6, &
math_Plain99to3333 math_Plain99to3333
use material, only: & use material, only: &
phasememberAt, &
phase_plasticity, & phase_plasticity, &
phase_plasticityInstance, &
material_phase, & material_phase, &
material_homog, & material_homog, &
temperature, & temperature, &
@ -490,7 +492,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
ho, & !< homogenization ho, & !< homogenization
tme !< thermal member position tme !< thermal member position
integer(pInt) :: & integer(pInt) :: &
i, j i, j, instance, of
ho = material_homog(ip,el) ho = material_homog(ip,el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
@ -509,8 +511,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) of = phasememberAt(ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
@ -818,6 +821,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_mul33x33, &
math_Mandel6to33, & math_Mandel6to33, &
math_Mandel33to6, & math_Mandel33to6, &
math_mul33x33 math_mul33x33
@ -825,6 +829,8 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: & use material, only: &
phasememberAt, &
phase_plasticityInstance, &
phase_plasticity, & phase_plasticity, &
phase_source, & phase_source, &
phase_Nsources, & phase_Nsources, &
@ -882,38 +888,41 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation) S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Mstar Mp
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
s !< counter in source loop s, & !< counter in source loop
instance, of
ho = material_homog( ip,el) ho = material_homog( ip,el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState (math_Mandel33to6(Mstar),ipc,ip,el) call plastic_isotropic_dotState (math_Mandel33to6(Mp),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(math_Mandel33to6(Mstar),ipc,ip,el) of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_phenopowerlaw_dotState(Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_dotState(math_Mandel33to6(Mstar),ipc,ip,el) call plastic_kinehardening_dotState(math_Mandel33to6(Mp),ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), & call plastic_dislotwin_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), & call plastic_disloucla_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (math_Mandel33to6(Mstar),FeArray,FpArray,temperature(ho)%p(tme), & call plastic_nonlocal_dotState (math_Mandel33to6(Mp),FeArray,FpArray,temperature(ho)%p(tme), &
subdt,subfracArray,ip,el) subdt,subfracArray,ip,el)
end select plasticityType end select plasticityType
@ -1027,13 +1036,18 @@ end subroutine constitutive_collectDeltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results !> @brief returns array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_postResults(S6, FeArray, ipc, ip, el) function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: &
math_Mandel6to33, &
math_mul33x33
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: & use material, only: &
phasememberAt, &
phase_plasticityInstance, &
plasticState, & plasticState, &
sourceState, & sourceState, &
phase_plasticity, & phase_plasticity, &
@ -1084,19 +1098,25 @@ function constitutive_postResults(S6, FeArray, ipc, ip, el)
real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults + & real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults + &
sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: & sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: &
constitutive_postResults constitutive_postResults
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
FeArray !< elastic deformation gradient FeArray !< elastic deformation gradient
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation) S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress
integer(pInt) :: & integer(pInt) :: &
startPos, endPos startPos, endPos
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
s !< counter in source loop s, of, instance !< counter in source loop
constitutive_postResults = 0.0_pReal constitutive_postResults = 0.0_pReal
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
ho = material_homog( ip,el) ho = material_homog( ip,el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
@ -1105,10 +1125,13 @@ function constitutive_postResults(S6, FeArray, ipc, ip, el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(S6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(S6,ipc,ip,el) plastic_isotropic_postResults(S6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_kinehardening_postResults(S6,ipc,ip,el) plastic_kinehardening_postResults(S6,ipc,ip,el)

4
src/crystallite.f90
View File

@ -3777,8 +3777,8 @@ function crystallite_postResults(ipc, ip, el)
c = c + 1_pInt c = c + 1_pInt
if (size(crystallite_postResults)-c > 0_pInt) & if (size(crystallite_postResults)-c > 0_pInt) &
crystallite_postResults(c+1:size(crystallite_postResults)) = & crystallite_postResults(c+1:size(crystallite_postResults)) = &
constitutive_postResults(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe, & constitutive_postResults(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc, ip, el) crystallite_Fe, ipc, ip, el)
end function crystallite_postResults end function crystallite_postResults

26
src/meshFEM.f90
View File

@ -8,10 +8,12 @@
!> results !> results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module mesh module mesh
#include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscis.h> #include <petsc/finclude/petscis.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use prec, only: pReal, pInt use prec, only: pReal, pInt
use PETScdmplex
use PETScdmda use PETScdmda
use PETScis use PETScis
@ -65,7 +67,7 @@ use PETScis
DM, public :: geomMesh DM, public :: geomMesh
integer(pInt), dimension(:), allocatable, public, protected :: & PetscInt, dimension(:), allocatable, public, protected :: &
mesh_boundaries mesh_boundaries
@ -88,15 +90,6 @@ use PETScis
mesh_FEM_build_ipCoordinates, & mesh_FEM_build_ipCoordinates, &
mesh_cellCenterCoordinates mesh_cellCenterCoordinates
external :: &
DMPlexCreateFromFile, &
DMPlexDistribute, &
DMPlexCopyCoordinates, &
DMGetStratumSize, &
DMPlexGetHeightStratum, &
DMPlexGetLabelValue, &
DMPlexSetLabelValue
contains contains
@ -151,14 +144,19 @@ subroutine mesh_init()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
! read in file
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr) call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
! get spatial dimension (2 or 3?)
call DMGetDimension(globalMesh,dimPlex,ierr) call DMGetDimension(globalMesh,dimPlex,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
write(6,*) 'dimension',dimPlex;flush(6)
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr) call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
! get number of IDs in face sets (for boundary conditions?)
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr) call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6)
call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
@ -174,7 +172,9 @@ subroutine mesh_init()
enddo enddo
if (nFaceSets > 0) call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr) if (nFaceSets > 0) call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr)
call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
! this read in function should ignore C and C++ style comments
! it is used for BC only?
if (worldrank == 0) then if (worldrank == 0) then
j = 0 j = 0
flag = .false. flag = .false.
@ -183,8 +183,8 @@ subroutine mesh_init()
read(FILEUNIT,'(a512)') line read(FILEUNIT,'(a512)') line
if (trim(line) == IO_EOF) exit ! skip empty lines if (trim(line) == IO_EOF) exit ! skip empty lines
if (trim(line) == '$Elements') then if (trim(line) == '$Elements') then
read(FILEUNIT,'(a512)') line read(FILEUNIT,'(a512)') line ! number of elements (ignore)
read(FILEUNIT,'(a512)') line read(FILEUNIT,'(a512)') line
flag = .true. flag = .true.
endif endif
if (trim(line) == '$EndElements') exit if (trim(line) == '$EndElements') exit

59
src/numerics.f90
View File

@ -127,12 +127,7 @@ module numerics
#ifdef FEM #ifdef FEM
integer(pInt), protected, public :: & integer(pInt), protected, public :: &
integrationOrder = 2_pInt, & !< order of quadrature rule required integrationOrder = 2_pInt, & !< order of quadrature rule required
structOrder = 2_pInt, & !< order of displacement shape functions structOrder = 2_pInt !< order of displacement shape functions
thermalOrder = 2_pInt, & !< order of temperature field shape functions
damageOrder = 2_pInt, & !< order of damage field shape functions
vacancyfluxOrder = 2_pInt, & !< order of vacancy concentration and chemical potential field shape functions
porosityOrder = 2_pInt, & !< order of porosity field shape functions
hydrogenfluxOrder = 2_pInt !< order of hydrogen concentration and chemical potential field shape functions
logical, protected, public :: & logical, protected, public :: &
BBarStabilisation = .false. BBarStabilisation = .false.
character(len=4096), protected, public :: & character(len=4096), protected, public :: &
@ -146,40 +141,7 @@ module numerics
&-mech_pc_type ml & &-mech_pc_type ml &
&-mech_mg_levels_ksp_type chebyshev & &-mech_mg_levels_ksp_type chebyshev &
&-mech_mg_levels_pc_type sor & &-mech_mg_levels_pc_type sor &
&-mech_pc_ml_nullspace user & &-mech_pc_ml_nullspace user ',&
&-damage_snes_type vinewtonrsls &
&-damage_snes_atol 1e-8 &
&-damage_ksp_type preonly &
&-damage_ksp_max_it 25 &
&-damage_pc_type cholesky &
&-damage_pc_factor_mat_solver_package mumps &
&-thermal_snes_type newtonls &
&-thermal_snes_linesearch_type cp &
&-thermal_ksp_type fgmres &
&-thermal_ksp_max_it 25 &
&-thermal_snes_atol 1e-3 &
&-thermal_pc_type hypre &
&-vacancy_snes_type newtonls &
&-vacancy_snes_linesearch_type cp &
&-vacancy_snes_atol 1e-9 &
&-vacancy_ksp_type fgmres &
&-vacancy_ksp_max_it 25 &
&-vacancy_pc_type ml &
&-vacancy_mg_levels_ksp_type chebyshev &
&-vacancy_mg_levels_pc_type sor &
&-porosity_snes_type newtonls &
&-porosity_snes_atol 1e-8 &
&-porosity_ksp_type fgmres &
&-porosity_ksp_max_it 25 &
&-porosity_pc_type hypre &
&-hydrogen_snes_type newtonls &
&-hydrogen_snes_linesearch_type cp &
&-hydrogen_snes_atol 1e-9 &
&-hydrogen_ksp_type fgmres &
&-hydrogen_ksp_max_it 25 &
&-hydrogen_pc_type ml &
&-hydrogen_mg_levels_ksp_type chebyshev &
&-hydrogen_mg_levels_pc_type sor ', &
petsc_options = '' petsc_options = ''
#endif #endif
@ -230,8 +192,6 @@ subroutine numerics_init
tag ,& tag ,&
line line
!$ character(len=6) DAMASK_NumThreadsString ! environment variable DAMASK_NUM_THREADS !$ character(len=6) DAMASK_NumThreadsString ! environment variable DAMASK_NUM_THREADS
external :: &
PETScErrorF ! is called in the CHKERRQ macro
#ifdef PETSc #ifdef PETSc
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr) call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
@ -458,16 +418,6 @@ subroutine numerics_init
integrationorder = IO_intValue(line,chunkPos,2_pInt) integrationorder = IO_intValue(line,chunkPos,2_pInt)
case ('structorder') case ('structorder')
structorder = IO_intValue(line,chunkPos,2_pInt) structorder = IO_intValue(line,chunkPos,2_pInt)
case ('thermalorder')
thermalorder = IO_intValue(line,chunkPos,2_pInt)
case ('damageorder')
damageorder = IO_intValue(line,chunkPos,2_pInt)
case ('vacancyfluxorder')
vacancyfluxOrder = IO_intValue(line,chunkPos,2_pInt)
case ('porosityorder')
porosityOrder = IO_intValue(line,chunkPos,2_pInt)
case ('hydrogenfluxorder')
hydrogenfluxOrder = IO_intValue(line,chunkPos,2_pInt)
case ('petsc_options') case ('petsc_options')
petsc_options = trim(line(chunkPos(4):)) petsc_options = trim(line(chunkPos(4):))
case ('bbarstabilisation') case ('bbarstabilisation')
@ -620,11 +570,6 @@ subroutine numerics_init
#ifdef FEM #ifdef FEM
write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder
write(6,'(a24,1x,i8)') ' structOrder: ',structOrder write(6,'(a24,1x,i8)') ' structOrder: ',structOrder
write(6,'(a24,1x,i8)') ' thermalOrder: ',thermalOrder
write(6,'(a24,1x,i8)') ' damageOrder: ',damageOrder
write(6,'(a24,1x,i8)') ' vacancyfluxOrder: ',vacancyfluxOrder
write(6,'(a24,1x,i8)') ' porosityOrder: ',porosityOrder
write(6,'(a24,1x,i8)') ' hydrogenfluxOrder: ',hydrogenfluxOrder
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_defaultOptions)//' '//trim(petsc_options) write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_defaultOptions)//' '//trim(petsc_options)
write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation
#endif #endif

822
src/plastic_phenopowerlaw.f90
View File

File diff suppressed because it is too large Load Diff

10
src/prec.f90
View File

@ -87,16 +87,6 @@ module prec
integer(pInt), pointer, dimension(:,:,:) :: p integer(pInt), pointer, dimension(:,:,:) :: p
end type end type
#ifdef FEM
type, public :: tOutputData
integer(pInt) :: &
sizeIpCells = 0_pInt , &
sizeResults = 0_pInt
real(pReal), allocatable, dimension(:,:) :: &
output !< output data
end type
#endif
public :: & public :: &
prec_init, & prec_init, &
dEq, & dEq, &

44
src/quit.f90 Normal file
View File

@ -0,0 +1,44 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief quit subroutine
!> @details exits the program and reports current time and duration. Exit code 0 signals
!> everything is fine. Exit code 1 signals an error, message according to IO_error. Exit code
!> 2 signals no severe problems, but some increments did not converge
!--------------------------------------------------------------------------------------------------
subroutine quit(stop_id)
#include <petsc/finclude/petscsys.h>
#ifdef _OPENMP
use MPI, only: &
MPI_finalize
#endif
use prec, only: &
pInt
use PetscSys
implicit none
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
integer(pInt) :: error = 0_pInt
PetscErrorCode :: ierr = 0
call PETScFinalize(ierr)
CHKERRQ(ierr)
#ifdef _OPENMP
call MPI_finalize(error)
if (error /= 0) write(6,'(a)') ' Error in MPI_finalize'
#endif
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
if (stop_id == 0_pInt .and. ierr == 0_pInt .and. error == 0_pInt) stop 0 ! normal termination
if (stop_id == 2_pInt .and. ierr == 0_pInt .and. error == 0_pInt) stop 2 ! not all incs converged
stop 1 ! error (message from IO_error)
end subroutine quit

12
src/spectral_damage.f90
View File

@ -50,8 +50,6 @@ module spectral_damage
spectral_damage_init, & spectral_damage_init, &
spectral_damage_solution, & spectral_damage_solution, &
spectral_damage_forward spectral_damage_forward
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains contains
@ -85,11 +83,6 @@ subroutine spectral_damage_init()
Vec :: uBound, lBound Vec :: uBound, lBound
PetscErrorCode :: ierr PetscErrorCode :: ierr
character(len=100) :: snes_type character(len=100) :: snes_type
external :: &
SNESSetOptionsPrefix, &
SNESGetType, &
DMDAGetCorners, &
DMDASNESSetFunctionLocal
write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>' write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press, ' write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press, '
@ -194,11 +187,6 @@ type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadC
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
external :: &
VecMin, &
VecMax, &
SNESSolve
spectral_damage_solution%converged =.false. spectral_damage_solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

10
src/spectral_mech_Basic.f90
View File

@ -66,8 +66,6 @@ module spectral_mech_basic
basic_init, & basic_init, &
basic_solution, & basic_solution, &
basic_forward basic_forward
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains contains
@ -119,11 +117,6 @@ subroutine basic_init
integer(pInt) :: proc integer(pInt) :: proc
character(len=1024) :: rankStr character(len=1024) :: rankStr
external :: &
SNESSetOptionsPrefix, &
SNESSetConvergenceTest, &
DMDASNESSetFunctionLocal
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:31–45, 2015' write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:31–45, 2015'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
@ -246,9 +239,6 @@ type(tSolutionState) function basic_solution(incInfoIn,timeinc,timeinc_old,stres
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
external :: &
SNESsolve
incInfo = incInfoIn incInfo = incInfoIn
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

10
src/spectral_mech_Polarisation.f90
View File

@ -73,8 +73,6 @@ module spectral_mech_Polarisation
Polarisation_init, & Polarisation_init, &
Polarisation_solution, & Polarisation_solution, &
Polarisation_forward Polarisation_forward
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains contains
@ -130,11 +128,6 @@ subroutine Polarisation_init
integer(pInt) :: proc integer(pInt) :: proc
character(len=1024) :: rankStr character(len=1024) :: rankStr
external :: &
SNESSetOptionsPrefix, &
SNESSetConvergenceTest, &
DMDASNESsetFunctionLocal
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:31–45, 2015' write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:31–45, 2015'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
@ -272,9 +265,6 @@ type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_ol
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
external :: &
SNESSolve
incInfo = incInfoIn incInfo = incInfoIn
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

12
src/spectral_thermal.f90
View File

@ -50,8 +50,6 @@ module spectral_thermal
spectral_thermal_init, & spectral_thermal_init, &
spectral_thermal_solution, & spectral_thermal_solution, &
spectral_thermal_forward spectral_thermal_forward
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains contains
@ -88,11 +86,6 @@ subroutine spectral_thermal_init
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr PetscErrorCode :: ierr
external :: &
SNESsetOptionsPrefix, &
DMDAgetCorners, &
DMDASNESsetFunctionLocal
write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>' write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press,' write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press,'
write(6,'(/,a)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018' write(6,'(/,a)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018'
@ -196,11 +189,6 @@ type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,load
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
external :: &
VecMin, &
VecMax, &
SNESSolve
spectral_thermal_solution%converged =.false. spectral_thermal_solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

4
src/spectral_utilities.f90
View File

@ -146,8 +146,6 @@ module spectral_utilities
FIELD_DAMAGE_ID FIELD_DAMAGE_ID
private :: & private :: &
utilities_getFreqDerivative utilities_getFreqDerivative
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains contains
@ -209,8 +207,6 @@ subroutine utilities_init()
scalarSize = 1_C_INTPTR_T, & scalarSize = 1_C_INTPTR_T, &
vecSize = 3_C_INTPTR_T, & vecSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T tensorSize = 9_C_INTPTR_T
external :: &
PetscOptionsInsertString
write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>' write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>'
write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity, 46:37–53, 2013' write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity, 46:37–53, 2013'