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

Merge branch '38-introduce-rudimentary-PETSc-based-FEM-solver' into 'development' 

Resolve "Introduce Rudimentary PETSc based FEM solver"

Closes #38

See merge request damask/DAMASK!41
This commit is contained in:
Martin Diehl 2018-10-01 21:10:03 +02:00
parent d5cf7d2144 14509876d4
commit 27bde05529
24 changed files with 324903 additions and 810 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

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

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")

387
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,12 +112,8 @@ 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(el = 1_pInt, ip = 1_pInt) call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
@ -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,133 +254,14 @@ 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
@ -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
@ -438,10 +293,9 @@ program DAMASK_FEM
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)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time
stepFraction = stepFraction + 1_pInt ! count step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop over sub incs ! report begin of new step
subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt)
time = time + timeinc ! forward time
stepFraction = stepFraction + 1_pInt
remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc
!--------------------------------------------------------------------------------------------------
! report begin of new increment
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 notConvergedCounter = notConvergedCounter + 1_pInt
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged' ' increment ', totalIncsCounter, ' NOT converged'
endif
notConvergedCounter = notConvergedCounter + 1_pInt
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
if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ...
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information
restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif endif
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ToDo first call to CPFEM_general will write? endif skipping
restartWrite = .true.
lastRestartWritten = inc
endif
else forwarding
time = time + timeinc
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)') &
convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', & notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/& real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!'
' %) increments converged!' flush(6)
endif call MPI_file_close(fileUnit,ierr)
if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged close(statUnit)
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_FEM 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(el = 1_pInt, ip = 1_pInt) call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
@ -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

81
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(:)
@ -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: &

236
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,24 +24,14 @@ 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, &
@ -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

22
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
@ -56,7 +58,7 @@ use PETScis
DM, public :: geomMesh DM, public :: geomMesh
integer(pInt), dimension(:), allocatable, public, protected :: & PetscInt, dimension(:), allocatable, public, protected :: &
mesh_boundaries mesh_boundaries
@ -89,15 +91,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
@ -154,14 +147,19 @@ subroutine mesh_init(ip,el)
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)
@ -178,6 +176,8 @@ subroutine mesh_init(ip,el)
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.
@ -186,7 +186,7 @@ subroutine mesh_init(ip,el)
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

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

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'