diff --git a/DAMASK_prerequisites.sh b/DAMASK_prerequisites.sh index d8f6824b9..12dd9bd07 100755 --- a/DAMASK_prerequisites.sh +++ b/DAMASK_prerequisites.sh @@ -9,14 +9,28 @@ OUTFILE="system_report.txt" echo =========================================== echo + Generating $OUTFILE echo + Send to damask@mpie.de for support +echo + view with \'cat $OUTFILE\' echo =========================================== - +function getDetails { +if which $1 &> /dev/null; then + echo ---------------------------------------------------------------------------------------------- + echo $1: + echo ---------------------------------------------------------------------------------------------- + echo + location: + which $1 + echo + $1 $2: + $1 $2 + echo -e '\n' +else + echo $ does not exist +fi +} # redirect STDOUT and STDERR to logfile # https://stackoverflow.com/questions/11229385/redirect-all-output-in-a-bash-script-when-using-set-x^ exec > $OUTFILE 2>&1 - + # directory, file is not a symlink by definition # https://stackoverflow.com/questions/59895/getting-the-source-directory-of-a-bash-script-from-within DAMASK_ROOT="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" @@ -28,13 +42,19 @@ echo echo ============================================================================================== echo DAMASK settings echo ============================================================================================== +echo ---------------------------------------------------------------------------------------------- echo DAMASK_ROOT: +echo ---------------------------------------------------------------------------------------------- echo $DAMASK_ROOT echo +echo ---------------------------------------------------------------------------------------------- echo Version: +echo ---------------------------------------------------------------------------------------------- cat VERSION echo +echo ---------------------------------------------------------------------------------------------- echo Settings in CONFIG: +echo ---------------------------------------------------------------------------------------------- cat CONFIG echo echo ============================================================================================== @@ -46,6 +66,8 @@ echo PATH: $PATH echo LD_LIBRARY_PATH: $LD_LIBRARY_PATH echo PYTHONPATH: $PYTHONPATH echo SHELL: $SHELL +echo PETSC_ARCH: $PETSC_ARCH +echo PETSC_DIR: $PETSC_DIR echo echo ============================================================================================== echo Python @@ -53,16 +75,14 @@ echo =========================================================================== DEFAULT_PYTHON=python2.7 for executable in python python2 python3 python2.7; do - if which $executable &> /dev/null; then - echo $executable version: $($executable --version 2>&1) - else - echo $executable does not exist - fi + getDetails $executable '--version' done -echo Details on $DEFAULT_PYTHON: $(ls -la $(which $DEFAULT_PYTHON)) -echo +echo ---------------------------------------------------------------------------------------------- +echo Details on $DEFAULT_PYTHON: +echo ---------------------------------------------------------------------------------------------- +echo $(ls -la $(which $DEFAULT_PYTHON)) for module in numpy scipy;do - echo ---------------------------------------------------------------------------------------------- + echo -e '\n----------------------------------------------------------------------------------------------' echo $module echo ---------------------------------------------------------------------------------------------- $DEFAULT_PYTHON -c "import $module; \ @@ -86,31 +106,36 @@ echo =========================================================================== echo GCC echo ============================================================================================== for executable in gcc g++ gfortran ;do - if which $executable &> /dev/null; then - echo $(which $executable) version: $($executable --version 2>&1) - else - echo $executable does not exist - fi + getDetails $executable '--version' done echo echo ============================================================================================== echo Intel Compiler Suite echo ============================================================================================== for executable in icc icpc ifort ;do - if which $executable &> /dev/null; then - echo $(which $executable) version: $($executable --version 2>&1) - else - echo $executable does not exist - fi + getDetails $executable '--version' done echo echo ============================================================================================== echo MPI Wrappers echo ============================================================================================== -for executable in mpicc mpiCC mpicxx mpicxx mpifort mpif90 mpif77; do - if which $executable &> /dev/null; then - echo $(which $executable) version: $($executable --show 2>&1) - else - echo $executable does not exist - fi +for executable in mpicc mpiCC mpic++ mpicpc mpicxx mpifort mpif90 mpif77; do + getDetails $executable '-show' done +echo +echo ============================================================================================== +echo MPI Launchers +echo ============================================================================================== +for executable in mpirun mpiexec; do + getDetails $executable '--version' +done +echo +echo ============================================================================================== +echo Abaqus +echo ============================================================================================== +cd installation/mods_Abaqus # to have the right environment file +for executable in abaqus abq2016 abq2017; do + getDetails $executable 'information=all' +done +cd ../.. + diff --git a/LICENSE b/LICENSE index f6371f259..630dc3a84 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright 2011-17 Max-Planck-Institut für Eisenforschung GmbH +Copyright 2011-18 Max-Planck-Institut für Eisenforschung GmbH DAMASK is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/VERSION b/VERSION index 35191bcab..becfecaca 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-1035-gd80a255 +v2.0.1-1107-g2c3eae6 diff --git a/examples/SpectralMethod/EshelbyInclusion/config/elastic_Ti.config b/examples/SpectralMethod/EshelbyInclusion/config/elastic_Ti.config new file mode 100644 index 000000000..b103eb112 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/config/elastic_Ti.config @@ -0,0 +1,7 @@ +elasticity hooke + +c11 160.0e9 +c12 90.0e9 +c13 66.0e9 +c33 181.7e9 +c44 46.5e9 diff --git a/examples/SpectralMethod/EshelbyInclusion/config/elastic_fullyAnisotropic.config b/examples/SpectralMethod/EshelbyInclusion/config/elastic_fullyAnisotropic.config new file mode 100644 index 000000000..700a37dbd --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/config/elastic_fullyAnisotropic.config @@ -0,0 +1,11 @@ +elasticity hooke + +c11 100.0e9 +c22 100.0e9 +c33 100.0e9 +c12 0.0e9 +c13 0.0e9 +c23 0.0e9 +c44 50.0e9 +c55 50.0e9 +c66 50.0e9 diff --git a/examples/SpectralMethod/EshelbyInclusion/config/elastic_isotropic.config b/examples/SpectralMethod/EshelbyInclusion/config/elastic_isotropic.config new file mode 100644 index 000000000..7574e9301 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/config/elastic_isotropic.config @@ -0,0 +1,4 @@ +elasticity hooke + +c11 100.0e9 +c12 0.0e9 diff --git a/examples/SpectralMethod/EshelbyInclusion/config/thermal.config b/examples/SpectralMethod/EshelbyInclusion/config/thermal.config new file mode 100644 index 000000000..50e50f3b3 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/config/thermal.config @@ -0,0 +1,3 @@ +reference_temperature 300.0 +specific_heat 1 +mass_density 1 diff --git a/examples/SpectralMethod/EshelbyInclusion/config/thermalExpansion_Ti.config b/examples/SpectralMethod/EshelbyInclusion/config/thermalExpansion_Ti.config new file mode 100644 index 000000000..5457be845 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/config/thermalExpansion_Ti.config @@ -0,0 +1,8 @@ +(kinematics) thermal_expansion +thermal_expansion11 9.5e-6 +thermal_expansion22 9.5e-6 +thermal_expansion33 5.6e-6 + +(source) thermal_externalheat +externalheat_time 0 500 500.001 1000 # 500 secs supplying 1 Watt, then removing 1 Watt +externalheat_rate 1 1 -1 -1 diff --git a/examples/SpectralMethod/EshelbyInclusion/config/thermalExpansion_fullyAnisotropic.config b/examples/SpectralMethod/EshelbyInclusion/config/thermalExpansion_fullyAnisotropic.config new file mode 100644 index 000000000..7190656d3 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/config/thermalExpansion_fullyAnisotropic.config @@ -0,0 +1,9 @@ +(kinematics) thermal_expansion +thermal_expansion11 5e-6 +thermal_expansion22 10e-6 +thermal_expansion33 20e-6 +lattice_structure orthorhombic + +(source) thermal_externalheat +externalheat_time 0 500 500.001 1000 # 500 secs supplying 1 Watt, then removing 1 Watt +externalheat_rate 1 1 -1 -1 diff --git a/examples/SpectralMethod/EshelbyInclusion/config/thermalExpansion_isotropic.config b/examples/SpectralMethod/EshelbyInclusion/config/thermalExpansion_isotropic.config new file mode 100644 index 000000000..5d039bec9 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/config/thermalExpansion_isotropic.config @@ -0,0 +1,6 @@ +(kinematics) thermal_expansion +thermal_expansion11 10e-6 + +(source) thermal_externalheat +externalheat_time 0 500 500.001 1000 # 500 secs supplying 1 Watt, then removing 1 Watt +externalheat_rate 1 1 -1 -1 diff --git a/examples/SpectralMethod/EshelbyInclusion/geom/Ti_Ti.geom b/examples/SpectralMethod/EshelbyInclusion/geom/Ti_Ti.geom new file mode 100644 index 000000000..28b49f368 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/geom/Ti_Ti.geom @@ -0,0 +1,427 @@ +9 header +geom_addPrimitive v2.0.1-1073-gc544fa1 -c 32 32 32 -d -16 -16 -16 inclusion.geom +geom_canvas v2.0.1-1073-gc544fa1 -o 16 16 16 -g 32 32 32 +geom_translate v2.0.1-1073-gc544fa1 -s 1,2,2,6 Ti_Ti.geom +geom_pack v2.0.1-1073-gc544fa1 Ti_Ti.geom isotropic_anisotropic.geom isotropic_isotropic.geom isotropic_rotated.geom +grid a 32 b 32 c 32 +size x 0.5 y 0.5 z 0.5 +origin x 0.25 y 0.25 z 0.25 +homogenization 1 +microstructures 2 +8623 of 2 +2 of 6 +29 of 2 +4 of 6 +27 of 2 +6 of 6 +26 of 2 +6 of 6 +27 of 2 +4 of 6 +29 of 2 +2 of 6 +798 of 2 +2 of 6 +28 of 2 +6 of 6 +25 of 2 +8 of 6 +24 of 2 +8 of 6 +23 of 2 +10 of 6 +22 of 2 +10 of 6 +23 of 2 +8 of 6 +24 of 2 +8 of 6 +25 of 2 +6 of 6 +28 of 2 +2 of 6 +701 of 2 +4 of 6 +26 of 2 +8 of 6 +23 of 2 +10 of 6 +22 of 2 +10 of 6 +21 of 2 +12 of 6 +20 of 2 +12 of 6 +20 of 2 +12 of 6 +20 of 2 +12 of 6 +21 of 2 +10 of 6 +22 of 2 +10 of 6 +23 of 2 +8 of 6 +26 of 2 +4 of 6 +637 of 2 +2 of 6 +27 of 2 +8 of 6 +23 of 2 +10 of 6 +21 of 2 +12 of 6 +20 of 2 +12 of 6 +20 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +20 of 2 +12 of 6 +20 of 2 +12 of 6 +21 of 2 +10 of 6 +23 of 2 +8 of 6 +27 of 2 +2 of 6 +604 of 2 +6 of 6 +24 of 2 +10 of 6 +21 of 2 +12 of 6 +20 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +20 of 2 +12 of 6 +21 of 2 +10 of 6 +24 of 2 +6 of 6 +572 of 2 +2 of 6 +27 of 2 +8 of 6 +23 of 2 +10 of 6 +21 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +17 of 2 +16 of 6 +16 of 2 +16 of 6 +17 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +21 of 2 +10 of 6 +23 of 2 +8 of 6 +27 of 2 +2 of 6 +541 of 2 +4 of 6 +26 of 2 +8 of 6 +22 of 2 +12 of 6 +20 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +17 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +17 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +20 of 2 +12 of 6 +22 of 2 +8 of 6 +26 of 2 +4 of 6 +539 of 2 +6 of 6 +24 of 2 +10 of 6 +21 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +17 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +17 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +21 of 2 +10 of 6 +24 of 2 +6 of 6 +538 of 2 +6 of 6 +24 of 2 +10 of 6 +21 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +17 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +17 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +21 of 2 +10 of 6 +24 of 2 +6 of 6 +539 of 2 +4 of 6 +26 of 2 +8 of 6 +22 of 2 +12 of 6 +20 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +17 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +16 of 2 +16 of 6 +17 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +20 of 2 +12 of 6 +22 of 2 +8 of 6 +26 of 2 +4 of 6 +541 of 2 +2 of 6 +27 of 2 +8 of 6 +23 of 2 +10 of 6 +21 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +17 of 2 +16 of 6 +16 of 2 +16 of 6 +17 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +21 of 2 +10 of 6 +23 of 2 +8 of 6 +27 of 2 +2 of 6 +572 of 2 +6 of 6 +24 of 2 +10 of 6 +21 of 2 +12 of 6 +20 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +20 of 2 +12 of 6 +21 of 2 +10 of 6 +24 of 2 +6 of 6 +604 of 2 +2 of 6 +27 of 2 +8 of 6 +23 of 2 +10 of 6 +21 of 2 +12 of 6 +20 of 2 +12 of 6 +20 of 2 +12 of 6 +19 of 2 +14 of 6 +18 of 2 +14 of 6 +19 of 2 +12 of 6 +20 of 2 +12 of 6 +20 of 2 +12 of 6 +21 of 2 +10 of 6 +23 of 2 +8 of 6 +27 of 2 +2 of 6 +637 of 2 +4 of 6 +26 of 2 +8 of 6 +23 of 2 +10 of 6 +22 of 2 +10 of 6 +21 of 2 +12 of 6 +20 of 2 +12 of 6 +20 of 2 +12 of 6 +20 of 2 +12 of 6 +21 of 2 +10 of 6 +22 of 2 +10 of 6 +23 of 2 +8 of 6 +26 of 2 +4 of 6 +701 of 2 +2 of 6 +28 of 2 +6 of 6 +25 of 2 +8 of 6 +24 of 2 +8 of 6 +23 of 2 +10 of 6 +22 of 2 +10 of 6 +23 of 2 +8 of 6 +24 of 2 +8 of 6 +25 of 2 +6 of 6 +28 of 2 +2 of 6 +798 of 2 +2 of 6 +29 of 2 +4 of 6 +27 of 2 +6 of 6 +26 of 2 +6 of 6 +27 of 2 +4 of 6 +29 of 2 +2 of 6 +8623 of 2 diff --git a/examples/SpectralMethod/EshelbyInclusion/geom/isotropic_anisotropic.geom b/examples/SpectralMethod/EshelbyInclusion/geom/isotropic_anisotropic.geom new file mode 100644 index 000000000..3d4d3eab5 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/geom/isotropic_anisotropic.geom @@ -0,0 +1,427 @@ +9 header +geom_addPrimitive v2.0.1-1073-gc544fa1 -c 32 32 32 -d -16 -16 -16 inclusion.geom +geom_canvas v2.0.1-1073-gc544fa1 -o 16 16 16 -g 32 32 32 +geom_translate v2.0.1-1073-gc544fa1 -s 2,4 isotropic_anisotropic.geom +geom_pack v2.0.1-1073-gc544fa1 Ti_Ti.geom isotropic_anisotropic.geom isotropic_isotropic.geom isotropic_rotated.geom +grid a 32 b 32 c 32 +size x 0.5 y 0.5 z 0.5 +origin x 0.25 y 0.25 z 0.25 +homogenization 1 +microstructures 2 +8623 of 1 +2 of 4 +29 of 1 +4 of 4 +27 of 1 +6 of 4 +26 of 1 +6 of 4 +27 of 1 +4 of 4 +29 of 1 +2 of 4 +798 of 1 +2 of 4 +28 of 1 +6 of 4 +25 of 1 +8 of 4 +24 of 1 +8 of 4 +23 of 1 +10 of 4 +22 of 1 +10 of 4 +23 of 1 +8 of 4 +24 of 1 +8 of 4 +25 of 1 +6 of 4 +28 of 1 +2 of 4 +701 of 1 +4 of 4 +26 of 1 +8 of 4 +23 of 1 +10 of 4 +22 of 1 +10 of 4 +21 of 1 +12 of 4 +20 of 1 +12 of 4 +20 of 1 +12 of 4 +20 of 1 +12 of 4 +21 of 1 +10 of 4 +22 of 1 +10 of 4 +23 of 1 +8 of 4 +26 of 1 +4 of 4 +637 of 1 +2 of 4 +27 of 1 +8 of 4 +23 of 1 +10 of 4 +21 of 1 +12 of 4 +20 of 1 +12 of 4 +20 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +20 of 1 +12 of 4 +20 of 1 +12 of 4 +21 of 1 +10 of 4 +23 of 1 +8 of 4 +27 of 1 +2 of 4 +604 of 1 +6 of 4 +24 of 1 +10 of 4 +21 of 1 +12 of 4 +20 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +20 of 1 +12 of 4 +21 of 1 +10 of 4 +24 of 1 +6 of 4 +572 of 1 +2 of 4 +27 of 1 +8 of 4 +23 of 1 +10 of 4 +21 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +17 of 1 +16 of 4 +16 of 1 +16 of 4 +17 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +21 of 1 +10 of 4 +23 of 1 +8 of 4 +27 of 1 +2 of 4 +541 of 1 +4 of 4 +26 of 1 +8 of 4 +22 of 1 +12 of 4 +20 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +17 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +17 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +20 of 1 +12 of 4 +22 of 1 +8 of 4 +26 of 1 +4 of 4 +539 of 1 +6 of 4 +24 of 1 +10 of 4 +21 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +17 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +17 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +21 of 1 +10 of 4 +24 of 1 +6 of 4 +538 of 1 +6 of 4 +24 of 1 +10 of 4 +21 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +17 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +17 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +21 of 1 +10 of 4 +24 of 1 +6 of 4 +539 of 1 +4 of 4 +26 of 1 +8 of 4 +22 of 1 +12 of 4 +20 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +17 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +16 of 1 +16 of 4 +17 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +20 of 1 +12 of 4 +22 of 1 +8 of 4 +26 of 1 +4 of 4 +541 of 1 +2 of 4 +27 of 1 +8 of 4 +23 of 1 +10 of 4 +21 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +17 of 1 +16 of 4 +16 of 1 +16 of 4 +17 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +21 of 1 +10 of 4 +23 of 1 +8 of 4 +27 of 1 +2 of 4 +572 of 1 +6 of 4 +24 of 1 +10 of 4 +21 of 1 +12 of 4 +20 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +20 of 1 +12 of 4 +21 of 1 +10 of 4 +24 of 1 +6 of 4 +604 of 1 +2 of 4 +27 of 1 +8 of 4 +23 of 1 +10 of 4 +21 of 1 +12 of 4 +20 of 1 +12 of 4 +20 of 1 +12 of 4 +19 of 1 +14 of 4 +18 of 1 +14 of 4 +19 of 1 +12 of 4 +20 of 1 +12 of 4 +20 of 1 +12 of 4 +21 of 1 +10 of 4 +23 of 1 +8 of 4 +27 of 1 +2 of 4 +637 of 1 +4 of 4 +26 of 1 +8 of 4 +23 of 1 +10 of 4 +22 of 1 +10 of 4 +21 of 1 +12 of 4 +20 of 1 +12 of 4 +20 of 1 +12 of 4 +20 of 1 +12 of 4 +21 of 1 +10 of 4 +22 of 1 +10 of 4 +23 of 1 +8 of 4 +26 of 1 +4 of 4 +701 of 1 +2 of 4 +28 of 1 +6 of 4 +25 of 1 +8 of 4 +24 of 1 +8 of 4 +23 of 1 +10 of 4 +22 of 1 +10 of 4 +23 of 1 +8 of 4 +24 of 1 +8 of 4 +25 of 1 +6 of 4 +28 of 1 +2 of 4 +798 of 1 +2 of 4 +29 of 1 +4 of 4 +27 of 1 +6 of 4 +26 of 1 +6 of 4 +27 of 1 +4 of 4 +29 of 1 +2 of 4 +8623 of 1 diff --git a/examples/SpectralMethod/EshelbyInclusion/geom/isotropic_isotropic.geom b/examples/SpectralMethod/EshelbyInclusion/geom/isotropic_isotropic.geom new file mode 100644 index 000000000..bf63f03d7 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/geom/isotropic_isotropic.geom @@ -0,0 +1,427 @@ +9 header +geom_addPrimitive v2.0.1-1073-gc544fa1 -c 32 32 32 -d -16 -16 -16 inclusion.geom +geom_canvas v2.0.1-1073-gc544fa1 -o 16 16 16 -g 32 32 32 +geom_translate v2.0.1-1073-gc544fa1 -s 2,3 isotropic_isotropic.geom +geom_pack v2.0.1-1073-gc544fa1 Ti_Ti.geom isotropic_anisotropic.geom isotropic_isotropic.geom isotropic_rotated.geom +grid a 32 b 32 c 32 +size x 0.5 y 0.5 z 0.5 +origin x 0.25 y 0.25 z 0.25 +homogenization 1 +microstructures 2 +8623 of 1 +2 of 3 +29 of 1 +4 of 3 +27 of 1 +6 of 3 +26 of 1 +6 of 3 +27 of 1 +4 of 3 +29 of 1 +2 of 3 +798 of 1 +2 of 3 +28 of 1 +6 of 3 +25 of 1 +8 of 3 +24 of 1 +8 of 3 +23 of 1 +10 of 3 +22 of 1 +10 of 3 +23 of 1 +8 of 3 +24 of 1 +8 of 3 +25 of 1 +6 of 3 +28 of 1 +2 of 3 +701 of 1 +4 of 3 +26 of 1 +8 of 3 +23 of 1 +10 of 3 +22 of 1 +10 of 3 +21 of 1 +12 of 3 +20 of 1 +12 of 3 +20 of 1 +12 of 3 +20 of 1 +12 of 3 +21 of 1 +10 of 3 +22 of 1 +10 of 3 +23 of 1 +8 of 3 +26 of 1 +4 of 3 +637 of 1 +2 of 3 +27 of 1 +8 of 3 +23 of 1 +10 of 3 +21 of 1 +12 of 3 +20 of 1 +12 of 3 +20 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +20 of 1 +12 of 3 +20 of 1 +12 of 3 +21 of 1 +10 of 3 +23 of 1 +8 of 3 +27 of 1 +2 of 3 +604 of 1 +6 of 3 +24 of 1 +10 of 3 +21 of 1 +12 of 3 +20 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +20 of 1 +12 of 3 +21 of 1 +10 of 3 +24 of 1 +6 of 3 +572 of 1 +2 of 3 +27 of 1 +8 of 3 +23 of 1 +10 of 3 +21 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +17 of 1 +16 of 3 +16 of 1 +16 of 3 +17 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +21 of 1 +10 of 3 +23 of 1 +8 of 3 +27 of 1 +2 of 3 +541 of 1 +4 of 3 +26 of 1 +8 of 3 +22 of 1 +12 of 3 +20 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +17 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +17 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +20 of 1 +12 of 3 +22 of 1 +8 of 3 +26 of 1 +4 of 3 +539 of 1 +6 of 3 +24 of 1 +10 of 3 +21 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +17 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +17 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +21 of 1 +10 of 3 +24 of 1 +6 of 3 +538 of 1 +6 of 3 +24 of 1 +10 of 3 +21 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +17 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +17 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +21 of 1 +10 of 3 +24 of 1 +6 of 3 +539 of 1 +4 of 3 +26 of 1 +8 of 3 +22 of 1 +12 of 3 +20 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +17 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +16 of 1 +16 of 3 +17 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +20 of 1 +12 of 3 +22 of 1 +8 of 3 +26 of 1 +4 of 3 +541 of 1 +2 of 3 +27 of 1 +8 of 3 +23 of 1 +10 of 3 +21 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +17 of 1 +16 of 3 +16 of 1 +16 of 3 +17 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +21 of 1 +10 of 3 +23 of 1 +8 of 3 +27 of 1 +2 of 3 +572 of 1 +6 of 3 +24 of 1 +10 of 3 +21 of 1 +12 of 3 +20 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +20 of 1 +12 of 3 +21 of 1 +10 of 3 +24 of 1 +6 of 3 +604 of 1 +2 of 3 +27 of 1 +8 of 3 +23 of 1 +10 of 3 +21 of 1 +12 of 3 +20 of 1 +12 of 3 +20 of 1 +12 of 3 +19 of 1 +14 of 3 +18 of 1 +14 of 3 +19 of 1 +12 of 3 +20 of 1 +12 of 3 +20 of 1 +12 of 3 +21 of 1 +10 of 3 +23 of 1 +8 of 3 +27 of 1 +2 of 3 +637 of 1 +4 of 3 +26 of 1 +8 of 3 +23 of 1 +10 of 3 +22 of 1 +10 of 3 +21 of 1 +12 of 3 +20 of 1 +12 of 3 +20 of 1 +12 of 3 +20 of 1 +12 of 3 +21 of 1 +10 of 3 +22 of 1 +10 of 3 +23 of 1 +8 of 3 +26 of 1 +4 of 3 +701 of 1 +2 of 3 +28 of 1 +6 of 3 +25 of 1 +8 of 3 +24 of 1 +8 of 3 +23 of 1 +10 of 3 +22 of 1 +10 of 3 +23 of 1 +8 of 3 +24 of 1 +8 of 3 +25 of 1 +6 of 3 +28 of 1 +2 of 3 +798 of 1 +2 of 3 +29 of 1 +4 of 3 +27 of 1 +6 of 3 +26 of 1 +6 of 3 +27 of 1 +4 of 3 +29 of 1 +2 of 3 +8623 of 1 diff --git a/examples/SpectralMethod/EshelbyInclusion/geom/isotropic_rotated.geom b/examples/SpectralMethod/EshelbyInclusion/geom/isotropic_rotated.geom new file mode 100644 index 000000000..66f8a8f1a --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/geom/isotropic_rotated.geom @@ -0,0 +1,427 @@ +9 header +geom_addPrimitive v2.0.1-1073-gc544fa1 -c 32 32 32 -d -16 -16 -16 inclusion.geom +geom_canvas v2.0.1-1073-gc544fa1 -o 16 16 16 -g 32 32 32 +geom_translate v2.0.1-1073-gc544fa1 -s 2,5 isotropic_rotated.geom +geom_pack v2.0.1-1073-gc544fa1 Ti_Ti.geom isotropic_anisotropic.geom isotropic_isotropic.geom isotropic_rotated.geom +grid a 32 b 32 c 32 +size x 0.5 y 0.5 z 0.5 +origin x 0.25 y 0.25 z 0.25 +homogenization 1 +microstructures 2 +8623 of 1 +2 of 5 +29 of 1 +4 of 5 +27 of 1 +6 of 5 +26 of 1 +6 of 5 +27 of 1 +4 of 5 +29 of 1 +2 of 5 +798 of 1 +2 of 5 +28 of 1 +6 of 5 +25 of 1 +8 of 5 +24 of 1 +8 of 5 +23 of 1 +10 of 5 +22 of 1 +10 of 5 +23 of 1 +8 of 5 +24 of 1 +8 of 5 +25 of 1 +6 of 5 +28 of 1 +2 of 5 +701 of 1 +4 of 5 +26 of 1 +8 of 5 +23 of 1 +10 of 5 +22 of 1 +10 of 5 +21 of 1 +12 of 5 +20 of 1 +12 of 5 +20 of 1 +12 of 5 +20 of 1 +12 of 5 +21 of 1 +10 of 5 +22 of 1 +10 of 5 +23 of 1 +8 of 5 +26 of 1 +4 of 5 +637 of 1 +2 of 5 +27 of 1 +8 of 5 +23 of 1 +10 of 5 +21 of 1 +12 of 5 +20 of 1 +12 of 5 +20 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +20 of 1 +12 of 5 +20 of 1 +12 of 5 +21 of 1 +10 of 5 +23 of 1 +8 of 5 +27 of 1 +2 of 5 +604 of 1 +6 of 5 +24 of 1 +10 of 5 +21 of 1 +12 of 5 +20 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +20 of 1 +12 of 5 +21 of 1 +10 of 5 +24 of 1 +6 of 5 +572 of 1 +2 of 5 +27 of 1 +8 of 5 +23 of 1 +10 of 5 +21 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +17 of 1 +16 of 5 +16 of 1 +16 of 5 +17 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +21 of 1 +10 of 5 +23 of 1 +8 of 5 +27 of 1 +2 of 5 +541 of 1 +4 of 5 +26 of 1 +8 of 5 +22 of 1 +12 of 5 +20 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +17 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +17 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +20 of 1 +12 of 5 +22 of 1 +8 of 5 +26 of 1 +4 of 5 +539 of 1 +6 of 5 +24 of 1 +10 of 5 +21 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +17 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +17 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +21 of 1 +10 of 5 +24 of 1 +6 of 5 +538 of 1 +6 of 5 +24 of 1 +10 of 5 +21 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +17 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +17 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +21 of 1 +10 of 5 +24 of 1 +6 of 5 +539 of 1 +4 of 5 +26 of 1 +8 of 5 +22 of 1 +12 of 5 +20 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +17 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +16 of 1 +16 of 5 +17 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +20 of 1 +12 of 5 +22 of 1 +8 of 5 +26 of 1 +4 of 5 +541 of 1 +2 of 5 +27 of 1 +8 of 5 +23 of 1 +10 of 5 +21 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +17 of 1 +16 of 5 +16 of 1 +16 of 5 +17 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +21 of 1 +10 of 5 +23 of 1 +8 of 5 +27 of 1 +2 of 5 +572 of 1 +6 of 5 +24 of 1 +10 of 5 +21 of 1 +12 of 5 +20 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +20 of 1 +12 of 5 +21 of 1 +10 of 5 +24 of 1 +6 of 5 +604 of 1 +2 of 5 +27 of 1 +8 of 5 +23 of 1 +10 of 5 +21 of 1 +12 of 5 +20 of 1 +12 of 5 +20 of 1 +12 of 5 +19 of 1 +14 of 5 +18 of 1 +14 of 5 +19 of 1 +12 of 5 +20 of 1 +12 of 5 +20 of 1 +12 of 5 +21 of 1 +10 of 5 +23 of 1 +8 of 5 +27 of 1 +2 of 5 +637 of 1 +4 of 5 +26 of 1 +8 of 5 +23 of 1 +10 of 5 +22 of 1 +10 of 5 +21 of 1 +12 of 5 +20 of 1 +12 of 5 +20 of 1 +12 of 5 +20 of 1 +12 of 5 +21 of 1 +10 of 5 +22 of 1 +10 of 5 +23 of 1 +8 of 5 +26 of 1 +4 of 5 +701 of 1 +2 of 5 +28 of 1 +6 of 5 +25 of 1 +8 of 5 +24 of 1 +8 of 5 +23 of 1 +10 of 5 +22 of 1 +10 of 5 +23 of 1 +8 of 5 +24 of 1 +8 of 5 +25 of 1 +6 of 5 +28 of 1 +2 of 5 +798 of 1 +2 of 5 +29 of 1 +4 of 5 +27 of 1 +6 of 5 +26 of 1 +6 of 5 +27 of 1 +4 of 5 +29 of 1 +2 of 5 +8623 of 1 diff --git a/examples/SpectralMethod/EshelbyInclusion/material.config b/examples/SpectralMethod/EshelbyInclusion/material.config new file mode 100644 index 000000000..83298680b --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/material.config @@ -0,0 +1,111 @@ +#-------------------# + +#-------------------# + +[direct] +type none # isostrain 1 grain + +thermal adiabatic # thermal strain (stress) induced mass transport +initialT 300.0 +(output) temperature + +#-------------------# + +#-------------------# + +[aLittleSomething] + +(output) texture +(output) f +(output) p +(output) fe +(output) fi +(output) fp + +#-------------------# + +#-------------------# + +#................. +[isotropic matrix] + +lattice_structure isotropic +plasticity none +{config/elastic_isotropic.config} +{config/thermal.config} + +#................. +[Ti matrix] + +lattice_structure hex +covera_ratio 1.587 +plasticity none +{config/elastic_Ti.config} +{config/thermal.config} + +#................. +[isotropic inclusion] + +lattice_structure isotropic +plasticity none +{config/elastic_isotropic.config} +{config/thermal.config} +{config/thermalExpansion_isotropic.config} + +#................. +[anisotropic inclusion] + +lattice_structure orthorhombic +plasticity none +{config/elastic_fullyAnisotropic.config} +{config/thermal.config} +{config/thermalExpansion_fullyAnisotropic.config} + +#................. +[Ti inclusion] + +lattice_structure hex +covera_ratio 1.587 +plasticity none +{config/elastic_Ti.config} +{config/thermal.config} +{config/thermalExpansion_Ti.config} + +#--------------------------# + +#--------------------------# + +[isotropic matrix] +crystallite 1 +(constituent) phase 1 texture 1 fraction 1.0 + +[Ti matrix] +crystallite 1 +(constituent) phase 2 texture 1 fraction 1.0 + +[isotropic inclusion] +crystallite 1 +(constituent) phase 3 texture 1 fraction 1.0 + +[anisotropic inclusion] +crystallite 1 +(constituent) phase 4 texture 1 fraction 1.0 + +[rotated inclusion] +crystallite 1 +(constituent) phase 4 texture 2 fraction 1.0 + +[Ti inclusion] +crystallite 1 +(constituent) phase 5 texture 1 fraction 1.0 + +#--------------------------# + +#--------------------------# + +[cube] +(gauss) phi1 0.0 Phi 0.0 phi2 0.0 scatter 0.0 fraction 1.0 + +[rotated] +(gauss) phi1 0.0 Phi 45.0 phi2 0.0 scatter 0.0 fraction 1.0 + diff --git a/examples/SpectralMethod/EshelbyInclusion/numerics.config b/examples/SpectralMethod/EshelbyInclusion/numerics.config new file mode 100644 index 000000000..191c204c8 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/numerics.config @@ -0,0 +1,5 @@ +#spectralsolver polarisation +spectralderivative fwbw_difference +err_div_tolrel 1e-3 +itmin 2 +petsc_options -mech_snes_type anderson -mech_snes_anderson_beta 1.0 -mech_snes_anderson_restart 10 -thermal_snes_type anderson -thermal_snes_anderson_beta 1.0 diff --git a/examples/SpectralMethod/EshelbyInclusion/runAll.sh b/examples/SpectralMethod/EshelbyInclusion/runAll.sh new file mode 100755 index 000000000..27630bc39 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/runAll.sh @@ -0,0 +1,56 @@ +#!/usr/local/bin/bash + +for geom in $(ls geom/*.geom) +do + base=${geom%.geom} + base=${base#geom/} + name=${base}_thermal + vtr=${base}.vtr + + [[ -f ${name}.spectralOut ]] || \ + DAMASK_spectral \ + --workingdir ./ \ + --load thermal.load \ + --geom $geom \ + > ${name}.out + + if [ ! -f postProc/${name}_inc10.txt ] + then + postResults ${name}.spectralOut \ + --ho temperature \ + --cr f,fe,fi,fp,p \ + --split \ + --separation x,y,z \ + + addCauchy postProc/${name}_inc*.txt \ + + addDeviator postProc/${name}_inc*.txt \ + --spherical \ + --tensor p,Cauchy \ + + addDisplacement postProc/${name}_inc*.txt \ + --nodal \ + + fi + + geom_check ${geom} + + for inc in {00..10} + do + echo "generating postProc/${name}_inc${inc}.vtr" + cp geom/${vtr} postProc/${name}_inc${inc}.vtr + vtk_addRectilinearGridData \ + postProc/${name}_inc${inc}.txt \ + --inplace \ + --vtk postProc/${name}_inc${inc}.vtr \ + --data 'sph(p)','sph(Cauchy)',temperature \ + --tensor f,fe,fi,fp,p,Cauchy \ + + vtk_addRectilinearGridData \ + postProc/${name}_inc${inc}_nodal.txt \ + --inplace \ + --vtk postProc/${name}_inc${inc}.vtr \ + --data 'avg(f).pos','fluct(f).pos' \ + + done +done diff --git a/examples/SpectralMethod/EshelbyInclusion/thermal.load b/examples/SpectralMethod/EshelbyInclusion/thermal.load new file mode 100644 index 000000000..a13aaa9e3 --- /dev/null +++ b/examples/SpectralMethod/EshelbyInclusion/thermal.load @@ -0,0 +1 @@ +Fdot 0 0 0 0 0 0 0 0 0 stress * * * * * * * * * time 1000 incs 10 diff --git a/examples/SpectralMethod/20grains.seeds b/examples/SpectralMethod/Polycrystal/20grains.seeds similarity index 100% rename from examples/SpectralMethod/20grains.seeds rename to examples/SpectralMethod/Polycrystal/20grains.seeds diff --git a/examples/SpectralMethod/20grains16x16x16.geom b/examples/SpectralMethod/Polycrystal/20grains16x16x16.geom similarity index 100% rename from examples/SpectralMethod/20grains16x16x16.geom rename to examples/SpectralMethod/Polycrystal/20grains16x16x16.geom diff --git a/examples/SpectralMethod/20grains32x32x32.geom b/examples/SpectralMethod/Polycrystal/20grains32x32x32.geom similarity index 100% rename from examples/SpectralMethod/20grains32x32x32.geom rename to examples/SpectralMethod/Polycrystal/20grains32x32x32.geom diff --git a/examples/SpectralMethod/20grains64x64x64.geom b/examples/SpectralMethod/Polycrystal/20grains64x64x64.geom similarity index 100% rename from examples/SpectralMethod/20grains64x64x64.geom rename to examples/SpectralMethod/Polycrystal/20grains64x64x64.geom diff --git a/examples/SpectralMethod/material.config b/examples/SpectralMethod/Polycrystal/material.config similarity index 100% rename from examples/SpectralMethod/material.config rename to examples/SpectralMethod/Polycrystal/material.config diff --git a/examples/SpectralMethod/numerics.config b/examples/SpectralMethod/Polycrystal/numerics.config similarity index 100% rename from examples/SpectralMethod/numerics.config rename to examples/SpectralMethod/Polycrystal/numerics.config diff --git a/examples/SpectralMethod/shearXY.load b/examples/SpectralMethod/Polycrystal/shearXY.load similarity index 100% rename from examples/SpectralMethod/shearXY.load rename to examples/SpectralMethod/Polycrystal/shearXY.load diff --git a/examples/SpectralMethod/shearZX.load b/examples/SpectralMethod/Polycrystal/shearZX.load similarity index 100% rename from examples/SpectralMethod/shearZX.load rename to examples/SpectralMethod/Polycrystal/shearZX.load diff --git a/examples/SpectralMethod/tensionX.load b/examples/SpectralMethod/Polycrystal/tensionX.load similarity index 100% rename from examples/SpectralMethod/tensionX.load rename to examples/SpectralMethod/Polycrystal/tensionX.load diff --git a/processing/post/addCalculation.py b/processing/post/addCalculation.py index 85c52621b..024e81cb6 100755 --- a/processing/post/addCalculation.py +++ b/processing/post/addCalculation.py @@ -2,7 +2,7 @@ # -*- coding: UTF-8 no BOM -*- import os,re,sys,collections -import math # noqa +import math,scipy,scipy.linalg # noqa import numpy as np from optparse import OptionParser import damask diff --git a/processing/pre/geom_vicinityOffset.py b/processing/pre/geom_vicinityOffset.py index 263b4e8d9..d4eb8e097 100755 --- a/processing/pre/geom_vicinityOffset.py +++ b/processing/pre/geom_vicinityOffset.py @@ -29,9 +29,14 @@ parser.add_option('-m', '--microstructureoffset', type = 'int', metavar = 'int', help = 'offset (positive or negative) for tagged microstructure indices. '+ '"0" selects maximum microstructure index [%default]') +parser.add_option('-n', '--nonperiodic', + dest = 'mode', + action = 'store_const', const = 'nearest', + help = 'assume geometry to be non-periodic') parser.set_defaults(vicinity = 1, offset = 0, + mode = 'wrap', ) (options, filenames) = parser.parse_args() @@ -79,8 +84,8 @@ for name in filenames: if options.offset == 0: options.offset = microstructure.max() - microstructure = np.where(ndimage.filters.maximum_filter(microstructure,size=1+2*options.vicinity,mode='wrap') == - ndimage.filters.minimum_filter(microstructure,size=1+2*options.vicinity,mode='wrap'), + microstructure = np.where(ndimage.filters.maximum_filter(microstructure,size=1+2*options.vicinity,mode=options.mode) == + ndimage.filters.minimum_filter(microstructure,size=1+2*options.vicinity,mode=options.mode), microstructure, microstructure + options.offset) newInfo['microstructures'] = microstructure.max() diff --git a/processing/pre/hybridIA_linODFsampling.py b/processing/pre/hybridIA_linODFsampling.py index 128d054e1..e4735674a 100755 --- a/processing/pre/hybridIA_linODFsampling.py +++ b/processing/pre/hybridIA_linODFsampling.py @@ -251,7 +251,7 @@ for name in filenames: table.head_read() errors = [] - labels = ['phi1','Phi','phi2','intensity'] + labels = ['1_euler','2_euler','3_euler','intensity'] for i,index in enumerate(table.label_index(labels)): if index < 0: errors.append('label {} not present.'.format(labels[i])) diff --git a/processing/pre/seeds_fromRandom.py b/processing/pre/seeds_fromRandom.py index 0df292ed2..77bacdc8e 100755 --- a/processing/pre/seeds_fromRandom.py +++ b/processing/pre/seeds_fromRandom.py @@ -166,28 +166,15 @@ for name in filenames: grainEuler[2,:] *= 360.0 # phi_2 is uniformly distributed if not options.selective: - seeds = np.array([]) - - while len(seeds) < options.N: - - theSeeds = np.zeros((options.N,3),dtype=float) # seed positions array - gridpoints = random.sample(range(gridSize),options.N) # choose first N from random permutation of grid positions - - theSeeds[:,0] = (np.mod(gridpoints ,options.grid[0])\ - +np.random.random(options.N)) /options.grid[0] - theSeeds[:,1] = (np.mod(gridpoints// options.grid[0] ,options.grid[1])\ - +np.random.random(options.N)) /options.grid[1] - theSeeds[:,2] = (np.mod(gridpoints//(options.grid[1]*options.grid[0]),options.grid[2])\ - +np.random.random(options.N)) /options.grid[2] - - goodSeeds = theSeeds[np.all(theSeeds<=options.fraction,axis=1)] # pick seeds within threshold fraction - seeds = goodSeeds if len(seeds) == 0 else np.vstack((seeds,goodSeeds)) - if len(seeds) > options.N: seeds = seeds[:min(options.N,len(seeds))] - - seeds = seeds.T # switch layout to point index as last index - + n = np.maximum(np.ones(3),np.array(options.grid*options.fraction), + dtype=int,casting='unsafe') # find max grid indices within fraction + meshgrid = np.meshgrid(*map(np.arange,n),indexing='ij') # create a meshgrid within fraction + coords = np.vstack((meshgrid[0],meshgrid[1],meshgrid[2])).reshape(3,n.prod()).T # assemble list of 3D coordinates + seeds = ((random.sample(coords,options.N)+np.random.random(options.N*3).reshape(options.N,3))\ + / \ + (n/options.fraction)).T # pick options.N of those, rattle position, + # and rescale to fall within fraction else: - seeds = np.zeros((options.N,3),dtype=float) # seed positions array seeds[0] = np.random.random(3)*options.grid/max(options.grid) i = 1 # start out with one given point diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index e34b5baa8..66aa11433 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -113,7 +113,7 @@ end subroutine CPFEM_initAll !> @brief allocate the arrays defined in module CPFEM and initialize them !-------------------------------------------------------------------------------------------------- subroutine CPFEM_init -#ifdef __GFORTRAN__ +#if __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index a16aee54f..a89bfc294 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -82,7 +82,7 @@ end subroutine CPFEM_initAll !> @brief allocate the arrays defined in module CPFEM and initialize them !-------------------------------------------------------------------------------------------------- subroutine CPFEM_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -140,36 +140,28 @@ subroutine CPFEM_init write(rankStr,'(a1,i0)')'_',worldrank call IO_read_intFile(777,'recordedPhase'//trim(rankStr),modelName,size(material_phase)) - read (777,rec=1) material_phase - close (777) + read (777,rec=1) material_phase; close (777) call IO_read_realFile(777,'convergedF'//trim(rankStr),modelName,size(crystallite_F0)) - read (777,rec=1) crystallite_F0 - close (777) + read (777,rec=1) crystallite_F0; close (777) call IO_read_realFile(777,'convergedFp'//trim(rankStr),modelName,size(crystallite_Fp0)) - read (777,rec=1) crystallite_Fp0 - close (777) + read (777,rec=1) crystallite_Fp0; close (777) call IO_read_realFile(777,'convergedFi'//trim(rankStr),modelName,size(crystallite_Fi0)) - read (777,rec=1) crystallite_Fi0 - close (777) + read (777,rec=1) crystallite_Fi0; close (777) call IO_read_realFile(777,'convergedLp'//trim(rankStr),modelName,size(crystallite_Lp0)) - read (777,rec=1) crystallite_Lp0 - close (777) + read (777,rec=1) crystallite_Lp0; close (777) call IO_read_realFile(777,'convergedLi'//trim(rankStr),modelName,size(crystallite_Li0)) - read (777,rec=1) crystallite_Li0 - close (777) + read (777,rec=1) crystallite_Li0; close (777) call IO_read_realFile(777,'convergeddPdF'//trim(rankStr),modelName,size(crystallite_dPdF0)) - read (777,rec=1) crystallite_dPdF0 - close (777) + read (777,rec=1) crystallite_dPdF0; close (777) call IO_read_realFile(777,'convergedTstar'//trim(rankStr),modelName,size(crystallite_Tstar0_v)) - read (777,rec=1) crystallite_Tstar0_v - close (777) + read (777,rec=1) crystallite_Tstar0_v; close (777) call IO_read_realFile(777,'convergedStateConst'//trim(rankStr),modelName) m = 0_pInt @@ -296,25 +288,25 @@ if (restartWrite) then write(rankStr,'(a1,i0)')'_',worldrank call IO_write_jobRealFile(777,'recordedPhase'//trim(rankStr),size(material_phase)) - write (777,rec=1) material_phase; close (777) + write (777,rec=1) material_phase; close (777) call IO_write_jobRealFile(777,'convergedF'//trim(rankStr),size(crystallite_F0)) - write (777,rec=1) crystallite_F0; close (777) + write (777,rec=1) crystallite_F0; close (777) call IO_write_jobRealFile(777,'convergedFp'//trim(rankStr),size(crystallite_Fp0)) - write (777,rec=1) crystallite_Fp0; close (777) + write (777,rec=1) crystallite_Fp0; close (777) call IO_write_jobRealFile(777,'convergedFi'//trim(rankStr),size(crystallite_Fi0)) - write (777,rec=1) crystallite_Fi0; close (777) + write (777,rec=1) crystallite_Fi0; close (777) call IO_write_jobRealFile(777,'convergedLp'//trim(rankStr),size(crystallite_Lp0)) - write (777,rec=1) crystallite_Lp0; close (777) + write (777,rec=1) crystallite_Lp0; close (777) call IO_write_jobRealFile(777,'convergedLi'//trim(rankStr),size(crystallite_Li0)) - write (777,rec=1) crystallite_Li0; close (777) + write (777,rec=1) crystallite_Li0; close (777) call IO_write_jobRealFile(777,'convergeddPdF'//trim(rankStr),size(crystallite_dPdF0)) - write (777,rec=1) crystallite_dPdF0; close (777) + write (777,rec=1) crystallite_dPdF0; close (777) call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v)) write (777,rec=1) crystallite_Tstar0_v; close (777) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index 5cda40249..068aebbc6 100644 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -7,7 +7,7 @@ !> results !-------------------------------------------------------------------------------------------------- program DAMASK_spectral -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -78,8 +78,7 @@ program DAMASK_spectral FIELD_UNDEFINED_ID, & FIELD_MECH_ID, & FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID, & - utilities_calcPlasticity + FIELD_DAMAGE_ID use spectral_mech_Basic use spectral_mech_AL use spectral_mech_Polarisation @@ -157,19 +156,6 @@ program DAMASK_spectral MPI_finalize, & MPI_allreduce, & PETScFinalize -!-------------------------------------------------------------------------------------------------- -! variables related to stop criterion for yielding - real(pReal) :: plasticWorkOld, plasticWorkNew, & ! plastic work - eqTotalStrainOld, eqTotalStrainNew, & ! total equivalent strain - eqPlasticStrainOld, eqPlasticStrainNew, & ! total equivalent plastic strain - eqStressOld, eqStressNew , & ! equivalent stress - yieldStopValue - real(pReal), dimension(3,3) :: yieldStress,yieldStressOld,yieldStressNew, & - plasticStrainOld, plasticStrainNew, plasticStrainRate - integer(pInt) :: yieldResUnit = 0_pInt - integer(pInt) :: stressstrainUnit = 0_pInt - character(len=13) :: stopFlag - logical :: yieldStop, yieldStopSatisfied !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) @@ -227,8 +213,6 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! reading the load case and assign values to the allocated data structure - yieldStop = .False. - yieldStopSatisfied = .False. rewind(FILEUNIT) do line = IO_read(FILEUNIT) @@ -303,30 +287,10 @@ program DAMASK_spectral temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) enddo loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) - case('totalstrain') - yieldStop = .True. - stopFlag = 'totalStrain' - yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) - case('plasticstrain') - yieldStop = .True. - stopFlag = 'plasticStrain' - yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) - case('plasticwork') - yieldStop = .True. - stopFlag = 'plasticWork' - yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) end select enddo; enddo close(FILEUNIT) - if(yieldStop) then ! initialize variables related to yield stop - yieldStressNew = 0.0_pReal - plasticStrainNew = 0.0_pReal - eqStressNew = 0.0_pReal - eqTotalStrainNew = 0.0_pReal - eqPlasticStrainNew = 0.0_pReal - plasticWorkNew = 0.0_pReal - endif !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase @@ -387,7 +351,7 @@ program DAMASK_spectral if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency write(6,'(2x,a,i5)') 'output frequency: ', & loadCases(currentLoadCase)%outputfrequency - write(6,'(2x,a,i5,/)') 'restart frequency: ', & + write(6,'(2x,a,i5,/)') 'restart frequency: ', & loadCases(currentLoadCase)%restartfrequency if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message enddo checkLoadcases @@ -429,23 +393,23 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! write header of output file if (worldrank == 0) then - if (.not. appendToOutFile) then ! after restart, append to existing results file + if (.not. appendToOutFile) then ! after restart, append to existing results file open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& '.spectralOut',form='UNFORMATTED',status='REPLACE') - write(resUnit) 'load:', trim(loadCaseFile) ! ... and write header + write(resUnit) 'load:', trim(loadCaseFile) ! ... and write header write(resUnit) 'workingdir:', trim(getSolverWorkingDirectoryName()) write(resUnit) 'geometry:', trim(geometryFile) write(resUnit) 'grid:', grid write(resUnit) 'size:', geomSize write(resUnit) 'materialpoint_sizeResults:', materialpoint_sizeResults write(resUnit) 'loadcases:', size(loadCases) - write(resUnit) 'frequencies:', loadCases%outputfrequency ! one entry per LoadCase - write(resUnit) 'times:', loadCases%time ! one entry per LoadCase + write(resUnit) 'frequencies:', loadCases%outputfrequency ! one entry per LoadCase + write(resUnit) 'times:', loadCases%time ! one entry per LoadCase write(resUnit) 'logscales:', loadCases%logscale - write(resUnit) 'increments:', loadCases%incs ! one entry per LoadCase - write(resUnit) 'startingIncrement:', restartInc - 1_pInt ! start with writing out the previous inc + write(resUnit) 'increments:', loadCases%incs ! one entry per LoadCase + write(resUnit) 'startingIncrement:', restartInc ! start with writing out the previous inc write(resUnit) 'eoh' - close(resUnit) ! end of header + close(resUnit) ! end of header open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& '.sta',form='FORMATTED',status='REPLACE') write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file @@ -516,16 +480,15 @@ program DAMASK_spectral endif else ! not-1st currentLoadCase of logarithmic scale 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))& - -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal))) + -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1_pInt ,pReal)/& + real(loadCases(currentLoadCase)%incs ,pReal))) endif endif + timeinc = timeinc * real(subStepFactor,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 - - skipping: if (totalIncsCounter < restartInc) then ! not yet at restart inc? + skipping: if (totalIncsCounter <= restartInc) then ! not yet at restart inc? time = time + timeinc ! just advance time, skip already performed calculation guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference else skipping @@ -542,17 +505,20 @@ program DAMASK_spectral ! report begin of new step write(6,'(/,a)') ' ###########################################################################' write(6,'(1x,a,es12.5'//& - ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& - ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//& + ',a,'//IO_intOut(inc) //',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& + ',a,'//IO_intOut(stepFraction) //',a,'//IO_intOut(subStepFactor**cutBackLevel)//& ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') & 'Time', time, & - 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& - '-', stepFraction, '/', subStepFactor**cutBackLevel,& + 's: Increment ', inc,'/',loadCases(currentLoadCase)%incs,& + '-', stepFraction,'/',subStepFactor**cutBackLevel,& ' of load case ', currentLoadCase,'/',size(loadCases) - write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//& - ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') & - 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& - '-',stepFraction, '/', subStepFactor**cutBackLevel + write(incInfo,& + '(a,'//IO_intOut(totalIncsCounter)//& + ',a,'//IO_intOut(sum(loadCases%incs))//& + ',a,'//IO_intOut(stepFraction)//& + ',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') & + 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& + '-', stepFraction,'/',subStepFactor**cutBackLevel flush(6) !-------------------------------------------------------------------------------------------------- @@ -629,7 +595,6 @@ program DAMASK_spectral stagIter = stagIter + 1_pInt stagIterate = stagIter < stagItMax & .and. all(solres(:)%converged) & - .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration enddo @@ -675,12 +640,10 @@ program DAMASK_spectral endif; flush(6) if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency - if (worldrank == 0) & - write(6,'(1/,a)') ' ... writing results to file ......................................' - flush(6) + write(6,'(1/,a)') ' ... writing results to file ......................................' + flush(6) call materialpoint_postResults() call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr) - if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek') do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & @@ -699,88 +662,22 @@ program DAMASK_spectral lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write? endif - else forwarding - time = time + timeinc - guess = .true. - endif forwarding - - yieldCheck: if(yieldStop) then ! check if it yields or satisfies the certain stop condition - yieldStressOld = yieldStressNew - plasticStrainOld = plasticStrainNew - eqStressOld = eqStressNew - eqTotalStrainOld = eqTotalStrainNew - eqPlasticStrainOld = eqPlasticStrainNew - plasticWorkOld = plasticWorkNew - - call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, & - eqPlasticStrainNew, plasticWorkNew, loadCases(currentLoadCase)%rotation) - - if (worldrank == 0) then ! output the stress-strain curve to file if yield stop criterion is used - if ((currentLoadCase == 1_pInt) .and. (inc == 1_pInt)) then - open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.stressstrain',form='FORMATTED',status='REPLACE') - write(stressstrainUnit,*) 0.0_pReal, 0.0_pReal - write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew - close(stressstrainUnit) - else - open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.stressstrain',form='FORMATTED', position='APPEND', status='OLD') - write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew - close(stressstrainUnit) - endif - endif - - if(stopFlag == 'totalStrain') then - if(eqTotalStrainNew > yieldStopValue) then - yieldStress = yieldStressOld * (eqTotalStrainNew - yieldStopValue)/(eqTotalStrainNew - eqTotalStrainOld) & ! linear interpolation of stress values - + yieldStressNew * (yieldStopValue - eqTotalStrainOld)/(eqTotalStrainNew - eqTotalStrainOld) - plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) ! calculate plastic strain rate - yieldStopSatisfied = .True. - endif - elseif(stopFlag == 'plasticStrain') then - if(eqPlasticStrainNew > yieldStopValue) then - yieldStress = yieldStressOld * (eqPlasticStrainNew - yieldStopValue)/(eqPlasticStrainNew - eqPlasticStrainOld) & - + yieldStressNew * (yieldStopValue - eqPlasticStrainOld)/(eqPlasticStrainNew - eqPlasticStrainOld) - plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) - yieldStopSatisfied = .True. - endif - elseif(stopFlag == 'plasticWork') then - if(plasticWorkNew > yieldStopValue) then - yieldStress = yieldStressOld * (plasticWorkNew - yieldStopValue)/(plasticWorkNew - plasticWorkOld) & - + yieldStressNew * (yieldStopValue - plasticWorkOld)/(plasticWorkNew - plasticWorkOld) - plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) - yieldStopSatisfied = .True. - endif - endif - endif yieldCheck - - if (yieldStopSatisfied) then ! when yield, write the yield stress and strain rate to file and quit the job - if (worldrank == 0) then - open(newunit=yieldResUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.yield',form='FORMATTED',status='REPLACE') - do i = 1_pInt,3_pInt - write(yieldResUnit,*) (yieldStress(i,j), j=1,3) - enddo - do i = 1_pInt,3_pInt - write(yieldResUnit,*) (plasticStrainRate(i,j), j=1,3) - enddo - close(yieldResUnit) - call quit(0_pInt) - endif - endif + endif skipping enddo incLooping + enddo loadCaseLooping !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation write(6,'(/,a)') ' ###########################################################################' - write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & - notConvergedCounter + convergedCounter, ' (', & - real(convergedCounter, pReal)/& - real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & - ' %) increments converged!' + write(6,'(1x,'//IO_intOut(convergedCounter)//',a,'//IO_intOut(notConvergedCounter + convergedCounter)//',a,f5.1,a)') & + convergedCounter, ' out of ', & + notConvergedCounter + convergedCounter, ' (', & + real(convergedCounter, pReal)/& + real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & + ' %) increments converged!' flush(6) call MPI_file_close(resUnit,ierr) close(statUnit) diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 index 5fdb2ebf8..3853cb37f 100644 --- a/src/FEsolving.f90 +++ b/src/FEsolving.f90 @@ -43,7 +43,7 @@ contains !> solver the information is provided by the interface module !-------------------------------------------------------------------------------------------------- subroutine FE_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -79,20 +79,22 @@ subroutine FE_init #include "compilation_info.f90" modelName = getSolverJobName() + +#if defined(Spectral) || defined(FEM) + #ifdef Spectral restartInc = spectralRestartInc - if(restartInc <= 0_pInt) then - call IO_warning(warning_ID=34_pInt) - restartInc = 1_pInt - endif - restartRead = restartInc > 1_pInt ! only read in if "true" restart requested -#elif defined FEM +#endif +#ifdef FEM restartInc = FEMRestartInc - if(restartInc <= 0_pInt) then +#endif + + if(restartInc < 0_pInt) then call IO_warning(warning_ID=34_pInt) - restartInc = 1_pInt + restartInc = 0_pInt endif - restartRead = restartInc > 1_pInt + restartRead = restartInc > 0_pInt ! only read in if "true" restart requested + #else call IO_open_inputFile(FILEUNIT,modelName) rewind(FILEUNIT) @@ -131,19 +133,19 @@ subroutine FE_init do read (FILEUNIT,'(a1024)',END=200) line chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'restart' .and. & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) == 'file' .and. & - IO_lc(IO_stringValue(line,chunkPos,3_pInt)) == 'job' .and. & - IO_lc(IO_stringValue(line,chunkPos,4_pInt)) == 'id' ) & + if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'restart' & + .and. IO_lc(IO_stringValue(line,chunkPos,2_pInt)) == 'file' & + .and. IO_lc(IO_stringValue(line,chunkPos,3_pInt)) == 'job' & + .and. IO_lc(IO_stringValue(line,chunkPos,4_pInt)) == 'id' ) & modelName = IO_StringValue(line,chunkPos,6_pInt) enddo -#else +#else ! QUESTION: is this meaningful for the spectral/FEM case? call IO_open_inputFile(FILEUNIT,modelName) rewind(FILEUNIT) do read (FILEUNIT,'(a1024)',END=200) line chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt))=='*heading') then + if (IO_lc(IO_stringValue(line,chunkPos,1_pInt))=='*heading') then read (FILEUNIT,'(a1024)',END=200) line chunkPos = IO_stringPos(line) modelName = IO_StringValue(line,chunkPos,1_pInt) diff --git a/src/IO.f90 b/src/IO.f90 index d877379c7..9e8033f73 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -82,7 +82,7 @@ contains !> @brief only outputs revision number !-------------------------------------------------------------------------------------------------- subroutine IO_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -1268,11 +1268,11 @@ integer(pInt) function IO_countNumericalDataLines(fileUnit) line = IO_read(fileUnit) chunkPos = IO_stringPos(line) tmp = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) - if (verify(trim(tmp) ,"0123456789")/=0) then ! found keyword + if (verify(trim(tmp),'0123456789') == 0) then ! numerical values + IO_countNumericalDataLines = IO_countNumericalDataLines + 1_pInt + else line = IO_read(fileUnit, .true.) ! reset IO_read exit - else - IO_countNumericalDataLines = IO_countNumericalDataLines + 1_pInt endif enddo backspace(fileUnit) diff --git a/src/compilation_info.f90 b/src/compilation_info.f90 index 3fc12f1ee..f0ca4d4cc 100644 --- a/src/compilation_info.f90 +++ b/src/compilation_info.f90 @@ -1,8 +1,7 @@ -#ifdef __GFORTRAN__ - write(6,*) 'Compiled with ', compiler_version() !not supported by and ifort <= 15 (and old gfortran) +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 + write(6,*) 'Compiled with ', compiler_version() write(6,*) 'With options ', compiler_options() -#endif -#ifdef __INTEL_COMPILER +#else write(6,'(a,i4.4,a,i8.8)') ' Compiled with Intel fortran version ', __INTEL_COMPILER,& ', build date ', __INTEL_COMPILER_BUILD_DATE #endif diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 877f0e8fa..a62245f99 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -37,7 +37,7 @@ contains !> @brief allocates arrays pointing to array of the various constitutive modules !-------------------------------------------------------------------------------------------------- subroutine constitutive_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -902,7 +902,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) integer(pLongInt) :: & - tick, tock, & + tick = 0_pLongInt, & + tock = 0_pLongInt, & tickrate, & maxticks integer(pInt) :: & diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 6dce1b5f8..12bf19871 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -55,14 +55,14 @@ module crystallite crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_partionedLi0,& !< intermediate velocity grad at start of homog inc crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - crystallite_P, & !< 1st Piola-Kirchhoff stress per grain - crystallite_subF !< def grad to be reached at end of crystallite inc + crystallite_P !< 1st Piola-Kirchhoff stress per grain real(pReal), dimension(:,:,:,:,:), allocatable, private :: & crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_subFp0,& !< plastic def grad at start of crystallite inc crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step) crystallite_subFi0,& !< intermediate def grad at start of crystallite inc + crystallite_subF, & !< def grad to be reached at end of crystallite inc crystallite_subF0, & !< def grad at start of crystallite inc crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc crystallite_subLi0,& !< intermediate velocity grad at start of crystallite inc @@ -137,7 +137,7 @@ contains !> @brief allocates and initialize per grain variables !-------------------------------------------------------------------------------------------------- subroutine crystallite_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -193,7 +193,7 @@ subroutine crystallite_init c, & !< counter in integration point component loop i, & !< counter in integration point loop e, & !< counter in element loop - o, & !< counter in output loop + o = 0_pInt, & !< counter in output loop r, & !< counter in crystallite loop cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points @@ -1250,14 +1250,16 @@ subroutine crystallite_integrateStateRK4() use numerics, only: & numerics_integrationMode use debug, only: & +#ifdef DEBUG + debug_e, & + debug_i, & + debug_g, & +#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & - debug_e, & - debug_i, & - debug_g, & debug_StateLoopDistribution use FEsolving, only: & FEsolving_execElem, & @@ -1544,14 +1546,16 @@ subroutine crystallite_integrateStateRKCK45() use, intrinsic :: & IEEE_arithmetic use debug, only: & +#ifdef DEBUG + debug_e, & + debug_i, & + debug_g, & +#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & - debug_e, & - debug_i, & - debug_g, & debug_StateLoopDistribution use numerics, only: & rTol_crystalliteState, & @@ -2047,14 +2051,16 @@ subroutine crystallite_integrateStateAdaptiveEuler() use, intrinsic :: & IEEE_arithmetic use debug, only: & +#ifdef DEBUG + debug_e, & + debug_i, & + debug_g, & +#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & - debug_e, & - debug_i, & - debug_g, & debug_StateLoopDistribution use numerics, only: & rTol_crystalliteState, & @@ -2407,14 +2413,16 @@ subroutine crystallite_integrateStateEuler() use, intrinsic :: & IEEE_arithmetic use debug, only: & +#ifdef DEBUG + debug_e, & + debug_i, & + debug_g, & +#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & - debug_e, & - debug_i, & - debug_g, & debug_StateLoopDistribution use numerics, only: & numerics_integrationMode, & @@ -2630,9 +2638,11 @@ subroutine crystallite_integrateStateFPI() use, intrinsic :: & IEEE_arithmetic use debug, only: & +#ifdef DEBUG debug_e, & debug_i, & debug_g, & +#endif debug_level,& debug_crystallite, & debug_levelBasic, & @@ -3094,14 +3104,16 @@ logical function crystallite_stateJump(ipc,ip,el) IEEE_arithmetic use prec, only: & dNeq0 +#ifdef DEBUG use debug, only: & + debug_e, & + debug_i, & + debug_g, & debug_level, & debug_crystallite, & debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i, & - debug_g + debug_levelSelective +#endif use material, only: & plasticState, & sourceState, & @@ -3231,9 +3243,11 @@ logical function crystallite_integrateStress(& debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & +#ifdef DEBUG debug_e, & debug_i, & debug_g, & +#endif debug_cumLpCalls, & debug_cumLpTicks, & debug_StressLoopLpDistribution, & @@ -3259,7 +3273,9 @@ logical function crystallite_integrateStress(& math_Plain33to9, & math_Plain9to33, & math_Plain99to3333 +#ifdef DEBUG use mesh, only: mesh_element +#endif implicit none integer(pInt), intent(in):: el, & ! element index @@ -3323,8 +3339,8 @@ logical function crystallite_integrateStress(& p, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update - integer(pLongInt) tick, & - tock, & + integer(pLongInt) :: tick = 0_pLongInt, & + tock = 0_pLongInt, & tickrate, & maxticks diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 2f3014937..59956e7d1 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -41,7 +41,7 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- subroutine damage_local_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 4750f5949..a1f0f0cd5 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -16,7 +16,7 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- subroutine damage_none_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index cd6ba8a5b..92ad183e1 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -46,7 +46,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine damage_nonlocal_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/debug.f90 b/src/debug.f90 index cbd6e659a..feb1c9fe4 100644 --- a/src/debug.f90 +++ b/src/debug.f90 @@ -102,7 +102,7 @@ contains !> @brief reads in parameters from debug.config and allocates arrays !-------------------------------------------------------------------------------------------------- subroutine debug_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 2f4124c2b..7dbea41d5 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -71,7 +71,7 @@ contains !> @brief module initialization !-------------------------------------------------------------------------------------------------- subroutine homogenization_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/homogenization_RGC.f90 b/src/homogenization_RGC.f90 index 84cb594db..611268393 100644 --- a/src/homogenization_RGC.f90 +++ b/src/homogenization_RGC.f90 @@ -72,7 +72,7 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- subroutine homogenization_RGC_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/homogenization_isostrain.f90 b/src/homogenization_isostrain.f90 index 055bfbb46..b569e3737 100644 --- a/src/homogenization_isostrain.f90 +++ b/src/homogenization_isostrain.f90 @@ -49,7 +49,7 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- subroutine homogenization_isostrain_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/homogenization_none.f90 b/src/homogenization_none.f90 index 75d8bcd3a..b2d2f52a7 100644 --- a/src/homogenization_none.f90 +++ b/src/homogenization_none.f90 @@ -18,7 +18,7 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- subroutine homogenization_none_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/hydrogenflux_cahnhilliard.f90 b/src/hydrogenflux_cahnhilliard.f90 index 89479a9c9..569be97dc 100644 --- a/src/hydrogenflux_cahnhilliard.f90 +++ b/src/hydrogenflux_cahnhilliard.f90 @@ -51,7 +51,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine hydrogenflux_cahnhilliard_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/hydrogenflux_isoconc.f90 b/src/hydrogenflux_isoconc.f90 index bef2a8437..d1b13aa76 100644 --- a/src/hydrogenflux_isoconc.f90 +++ b/src/hydrogenflux_isoconc.f90 @@ -16,7 +16,7 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- subroutine hydrogenflux_isoconc_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index fffa26165..74af0a52d 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -51,7 +51,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine kinematics_cleavage_opening_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/kinematics_hydrogen_strain.f90 b/src/kinematics_hydrogen_strain.f90 index c3af7e2a2..f3ea4df38 100644 --- a/src/kinematics_hydrogen_strain.f90 +++ b/src/kinematics_hydrogen_strain.f90 @@ -41,7 +41,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine kinematics_hydrogen_strain_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 07b98aa23..ba38ac05b 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -51,7 +51,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine kinematics_slipplane_opening_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index e7cbca673..0de483d70 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -41,7 +41,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine kinematics_thermal_expansion_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -174,8 +174,12 @@ pure function kinematics_thermal_expansion_initialStrain(ipc, ip, el) offset = thermalMapping(homog)%p(ip,el) kinematics_thermal_expansion_initialStrain = & - (temperature(homog)%p(offset) - lattice_referenceTemperature(phase)) * & - lattice_thermalExpansion33(1:3,1:3,phase) + (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * & + lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient + (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * & + lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient + (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * & + lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient end function kinematics_thermal_expansion_initialStrain @@ -215,9 +219,16 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, TDot = temperatureRate(homog)%p(offset) TRef = lattice_referenceTemperature(phase) - Li = TDot* & - lattice_thermalExpansion33(1:3,1:3,phase)/ & - (1.0_pReal + lattice_thermalExpansion33(1:3,1:3,phase)*(T - TRef)) + Li = TDot * ( & + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient + + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient + + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient + ) / & + (1.0_pReal \ + + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. & + + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & + ) dLi_dTstar3333 = 0.0_pReal end subroutine kinematics_thermal_expansion_LiAndItsTangent diff --git a/src/kinematics_vacancy_strain.f90 b/src/kinematics_vacancy_strain.f90 index 9558f506d..227a86e0c 100644 --- a/src/kinematics_vacancy_strain.f90 +++ b/src/kinematics_vacancy_strain.f90 @@ -41,7 +41,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine kinematics_vacancy_strain_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/lattice.f90 b/src/lattice.f90 index 9635643e8..37393b82e 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1082,9 +1082,10 @@ module lattice lattice_nu, & lattice_trans_mu, & lattice_trans_nu + real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & ! with higher-order parameters (e.g. temperature-dependent) + lattice_thermalExpansion33 real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_thermalConductivity33, & - lattice_thermalExpansion33, & lattice_damageDiffusion33, & lattice_vacancyfluxDiffusion33, & lattice_vacancyfluxMobility33, & @@ -1243,7 +1244,7 @@ contains !> @brief Module initialization !-------------------------------------------------------------------------------------------------- subroutine lattice_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -1380,8 +1381,8 @@ subroutine lattice_init allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal) allocate(lattice_trans_C66(6,6,Nphases), source=0.0_pReal) allocate(lattice_trans_C3333(3,3,3,3,Nphases), source=0.0_pReal) + allocate(lattice_thermalExpansion33 (3,3,3,Nphases), source=0.0_pReal) ! constant, linear, quadratic coefficients allocate(lattice_thermalConductivity33 (3,3,Nphases), source=0.0_pReal) - allocate(lattice_thermalExpansion33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_damageDiffusion33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_vacancyfluxDiffusion33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_vacancyfluxMobility33 (3,3,Nphases), source=0.0_pReal) @@ -1545,11 +1546,17 @@ subroutine lattice_init case ('thermal_conductivity33') lattice_thermalConductivity33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) case ('thermal_expansion11') - lattice_thermalExpansion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) + do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T) + lattice_thermalExpansion33(1,1,i-1_pInt,section) = IO_floatValue(line,chunkPos,i) + enddo case ('thermal_expansion22') - lattice_thermalExpansion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) + do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T) + lattice_thermalExpansion33(2,2,i-1_pInt,section) = IO_floatValue(line,chunkPos,i) + enddo case ('thermal_expansion33') - lattice_thermalExpansion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) + do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T) + lattice_thermalExpansion33(3,3,i-1_pInt,section) = IO_floatValue(line,chunkPos,i) + enddo case ('specific_heat') lattice_specificHeat(section) = IO_floatValue(line,chunkPos,2_pInt) case ('vacancyformationenergy') @@ -1756,22 +1763,24 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) end select end select - lattice_thermalConductivity33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_thermalConductivity33(1:3,1:3,myPhase)) - lattice_thermalExpansion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_thermalExpansion33(1:3,1:3,myPhase)) - lattice_DamageDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_DamageDiffusion33(1:3,1:3,myPhase)) - lattice_vacancyfluxDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_vacancyfluxDiffusion33(1:3,1:3,myPhase)) - lattice_vacancyfluxMobility33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_vacancyfluxMobility33(1:3,1:3,myPhase)) - lattice_PorosityDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_PorosityDiffusion33(1:3,1:3,myPhase)) + forall (i = 1_pInt:3_pInt) & + lattice_thermalExpansion33 (1:3,1:3,i,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_thermalExpansion33 (1:3,1:3,i,myPhase)) + + lattice_thermalConductivity33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_thermalConductivity33 (1:3,1:3,myPhase)) + lattice_DamageDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_DamageDiffusion33 (1:3,1:3,myPhase)) + lattice_vacancyfluxDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_vacancyfluxDiffusion33 (1:3,1:3,myPhase)) + lattice_vacancyfluxMobility33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_vacancyfluxMobility33 (1:3,1:3,myPhase)) + lattice_PorosityDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_PorosityDiffusion33 (1:3,1:3,myPhase)) lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase)) - lattice_hydrogenfluxMobility33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_hydrogenfluxMobility33(1:3,1:3,myPhase)) + lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase)) + lattice_hydrogenfluxMobility33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_hydrogenfluxMobility33 (1:3,1:3,myPhase)) select case(lattice_structure(myPhase)) !-------------------------------------------------------------------------------------------------- diff --git a/src/material.f90 b/src/material.f90 index 551b9d981..25d115520 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -364,7 +364,7 @@ contains !> material.config !-------------------------------------------------------------------------------------------------- subroutine material_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/math.f90 b/src/math.f90 index af125e101..3c7f2d91c 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -174,7 +174,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine math_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -1743,13 +1743,12 @@ real(pReal) pure function math_EulerMisorientation(EulerA,EulerB) implicit none real(pReal), dimension(3), intent(in) :: EulerA,EulerB - real(pReal), dimension(3,3) :: r - real(pReal) :: tr + real(pReal) :: cosTheta - r = math_mul33x33(math_EulerToR(EulerB),transpose(math_EulerToR(EulerA))) + cosTheta = (math_trace33(math_mul33x33(math_EulerToR(EulerB), & + transpose(math_EulerToR(EulerA)))) - 1.0_pReal) * 0.5_pReal - tr = (math_trace33(r)-1.0_pReal)*0.4999999_pReal - math_EulerMisorientation = abs(0.5_pReal*PI-asin(tr)) + math_EulerMisorientation = acos(math_limit(cosTheta,-1.0_pReal,1.0_pReal)) end function math_EulerMisorientation @@ -1782,30 +1781,28 @@ function math_sampleGaussOri(center,noise) real(pReal), dimension(3), intent(in) :: center real(pReal) :: cosScatter,scatter real(pReal), dimension(3) :: math_sampleGaussOri, disturb - real(pReal), dimension(3), parameter :: ORIGIN = [0.0_pReal,0.0_pReal,0.0_pReal] + real(pReal), dimension(3), parameter :: ORIGIN = 0.0_pReal real(pReal), dimension(5) :: rnd - integer(pInt) :: i - if (abs(noise) < tol_math_check) then + noScatter: if (abs(noise) < tol_math_check) then math_sampleGaussOri = center - return - endif + else noScatter + ! Helming uses different distribution with Bessel functions + ! therefore the gauss scatter width has to be scaled differently + scatter = 0.95_pReal * noise + cosScatter = cos(scatter) -! Helming uses different distribution with Bessel functions -! therefore the gauss scatter width has to be scaled differently - scatter = 0.95_pReal * noise - cosScatter = cos(scatter) + do + call halton(5_pInt,rnd) + rnd(1:3) = 2.0_pReal*rnd(1:3)-1.0_pReal ! expand 1:3 to range [-1,+1] + disturb = [ scatter * rnd(1), & ! phi1 + sign(1.0_pReal,rnd(2))*acos(cosScatter+(1.0_pReal-cosScatter)*rnd(4)), & ! Phi + scatter * rnd(3)] ! phi2 + if (rnd(5) <= exp(-1.0_pReal*(math_EulerMisorientation(ORIGIN,disturb)/scatter)**2_pReal)) exit + enddo - do - call halton(5_pInt,rnd) - forall (i=1_pInt:3_pInt) rnd(i) = 2.0_pReal*rnd(i)-1.0_pReal ! expand 1:3 to range [-1,+1] - disturb = [ scatter * rnd(1), & ! phi1 - sign(1.0_pReal,rnd(2))*acos(cosScatter+(1.0_pReal-cosScatter)*rnd(4)), & ! Phi - scatter * rnd(2)] ! phi2 - if (rnd(5) <= exp(-1.0_pReal*(math_EulerMisorientation(ORIGIN,disturb)/scatter)**2_pReal)) exit - enddo - - math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center))) + math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center))) + endif noScatter end function math_sampleGaussOri @@ -2759,8 +2756,7 @@ pure function math_rotate_forward33(tensor,rot_tensor) real(pReal), dimension(3,3) :: math_rotate_forward33 real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor - math_rotate_forward33 = math_mul33x33(rot_tensor,& - math_mul33x33(tensor,math_transpose33(rot_tensor))) + math_rotate_forward33 = math_mul33x33(rot_tensor,math_mul33x33(tensor,transpose(rot_tensor))) end function math_rotate_forward33 @@ -2774,8 +2770,7 @@ pure function math_rotate_backward33(tensor,rot_tensor) real(pReal), dimension(3,3) :: math_rotate_backward33 real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor - math_rotate_backward33 = math_mul33x33(math_transpose33(rot_tensor),& - math_mul33x33(tensor,rot_tensor)) + math_rotate_backward33 = math_mul33x33(transpose(rot_tensor),math_mul33x33(tensor,rot_tensor)) end function math_rotate_backward33 @@ -2796,8 +2791,8 @@ pure function math_rotate_forward3333(tensor,rot_tensor) do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt; do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt; do o = 1_pInt,3_pInt; do p = 1_pInt,3_pInt math_rotate_forward3333(i,j,k,l) = math_rotate_forward3333(i,j,k,l) & - + rot_tensor(m,i) * rot_tensor(n,j) & - * rot_tensor(o,k) * rot_tensor(p,l) * tensor(m,n,o,p) + + rot_tensor(i,m) * rot_tensor(j,n) & + * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p) enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo end function math_rotate_forward3333 diff --git a/src/mesh.f90 b/src/mesh.f90 index d7d0f8c06..74bb55a3b 100644 --- a/src/mesh.f90 +++ b/src/mesh.f90 @@ -379,30 +379,30 @@ module mesh ],pInt) - integer(pInt), dimension(FE_Nelemtypes), parameter, private :: MESH_VTKELEMTYPE = & - int([ & - 5, & ! element 6 (2D 3node 1ip) - 22, & ! element 125 (2D 6node 3ip) - 9, & ! element 11 (2D 4node 4ip) - 23, & ! element 27 (2D 8node 9ip) - 23, & ! element 54 (2D 8node 4ip) - 10, & ! element 134 (3D 4node 1ip) - 10, & ! element 157 (3D 5node 4ip) - 24, & ! element 127 (3D 10node 4ip) - 13, & ! element 136 (3D 6node 6ip) - 12, & ! element 117 (3D 8node 1ip) - 12, & ! element 7 (3D 8node 8ip) - 25, & ! element 57 (3D 20node 8ip) - 25 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer(pInt), dimension(FE_Ncelltypes), parameter, private :: MESH_VTKCELLTYPE = & - int([ & - 5, & ! (2D 3node) - 9, & ! (2D 4node) - 10, & ! (3D 4node) - 12 & ! (3D 8node) - ],pInt) +! integer(pInt), dimension(FE_Nelemtypes), parameter, private :: MESH_VTKELEMTYPE = & +! int([ & +! 5, & ! element 6 (2D 3node 1ip) +! 22, & ! element 125 (2D 6node 3ip) +! 9, & ! element 11 (2D 4node 4ip) +! 23, & ! element 27 (2D 8node 9ip) +! 23, & ! element 54 (2D 8node 4ip) +! 10, & ! element 134 (3D 4node 1ip) +! 10, & ! element 157 (3D 5node 4ip) +! 24, & ! element 127 (3D 10node 4ip) +! 13, & ! element 136 (3D 6node 6ip) +! 12, & ! element 117 (3D 8node 1ip) +! 12, & ! element 7 (3D 8node 8ip) +! 25, & ! element 57 (3D 20node 8ip) +! 25 & ! element 21 (3D 20node 27ip) +! ],pInt) +! +! integer(pInt), dimension(FE_Ncelltypes), parameter, private :: MESH_VTKCELLTYPE = & +! int([ & +! 5, & ! (2D 3node) +! 9, & ! (2D 4node) +! 10, & ! (3D 4node) +! 12 & ! (3D 8node) +! ],pInt) public :: & @@ -477,7 +477,7 @@ contains !! Order and routines strongly depend on type of solver !-------------------------------------------------------------------------------------------------- subroutine mesh_init(ip,el) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -2848,16 +2848,18 @@ use IO, only: & implicit none integer(pInt), intent(in) :: fileUnit -#ifndef Spectral +#ifdef Spectral + mesh_periodicSurface = .true. + + end subroutine mesh_get_damaskOptions + +#else + integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) chunk, Nchunks character(len=300) :: line, damaskOption, v character(len=300) :: keyword -#endif -#ifdef Spectral - mesh_periodicSurface = .true. -#else mesh_periodicSurface = .false. #ifdef Marc4DAMASK keyword = '$damask' @@ -2886,9 +2888,9 @@ use IO, only: & enddo 610 FORMAT(A300) -#endif 620 end subroutine mesh_get_damaskOptions +#endif !-------------------------------------------------------------------------------------------------- diff --git a/src/numerics.f90 b/src/numerics.f90 index 8392ac61c..c854d9d2b 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -198,7 +198,7 @@ contains ! a sanity check !-------------------------------------------------------------------------------------------------- subroutine numerics_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index c02a7c4d4..9c0a6c494 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -119,7 +119,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_disloUCLA_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index d6c73e8f3..c7aaf5400 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -198,7 +198,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index f1ad909ba..4a9d70061 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -90,7 +90,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index 839a4fa9f..140754556 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -26,7 +26,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_none_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index facc76ff1..64520fdc6 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -123,7 +123,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_phenopowerlaw_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/porosity_none.f90 b/src/porosity_none.f90 index 2bca99384..c273baf3b 100644 --- a/src/porosity_none.f90 +++ b/src/porosity_none.f90 @@ -16,7 +16,7 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- subroutine porosity_none_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/porosity_phasefield.f90 b/src/porosity_phasefield.f90 index 3f1c853a4..6ab7263e4 100644 --- a/src/porosity_phasefield.f90 +++ b/src/porosity_phasefield.f90 @@ -48,7 +48,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine porosity_phasefield_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/prec.f90 b/src/prec.f90 index 912a02533..f35735780 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -111,7 +111,7 @@ contains !> @brief reporting precision !-------------------------------------------------------------------------------------------------- subroutine prec_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index cad6bf1e4..22236a636 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -63,7 +63,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoBrittle_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 959e62e26..944a65918 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -67,7 +67,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoDuctile_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index bd1026765..b9fb2c22c 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -53,7 +53,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoBrittle_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 5d7e4f862..ed08e0a41 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -53,7 +53,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoDuctile_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 99c41f062..7a4e85c75 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -39,7 +39,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_thermal_dissipation_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index 6b015689a..2907ddf85 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -45,7 +45,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/source_vacancy_irradiation.f90 b/src/source_vacancy_irradiation.f90 index 8f24b39be..566d97e68 100644 --- a/src/source_vacancy_irradiation.f90 +++ b/src/source_vacancy_irradiation.f90 @@ -41,7 +41,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_vacancy_irradiation_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/source_vacancy_phenoplasticity.f90 b/src/source_vacancy_phenoplasticity.f90 index 26c3ae828..8834a067a 100644 --- a/src/source_vacancy_phenoplasticity.f90 +++ b/src/source_vacancy_phenoplasticity.f90 @@ -39,7 +39,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_vacancy_phenoplasticity_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/source_vacancy_thermalfluc.f90 b/src/source_vacancy_thermalfluc.f90 index e5d3b0574..91047fbf2 100644 --- a/src/source_vacancy_thermalfluc.f90 +++ b/src/source_vacancy_thermalfluc.f90 @@ -41,7 +41,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_vacancy_thermalfluc_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/spectral_damage.f90 b/src/spectral_damage.f90 index 727659870..1ac3c4c73 100644 --- a/src/spectral_damage.f90 +++ b/src/spectral_damage.f90 @@ -60,7 +60,7 @@ contains !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine spectral_damage_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/spectral_interface.f90 b/src/spectral_interface.f90 index 3c8489d04..b45e8316c 100644 --- a/src/spectral_interface.f90 +++ b/src/spectral_interface.f90 @@ -15,7 +15,7 @@ module DAMASK_interface private #include logical, public, protected :: appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding) - integer(pInt), public, protected :: spectralRestartInc = 1_pInt !< Increment at which calculation starts + integer(pInt), public, protected :: spectralRestartInc = 0_pInt !< Increment at which calculation starts character(len=1024), public, protected :: & geometryFile = '', & !< parameter given for geometry file loadCaseFile = '' !< parameter given for load case file @@ -83,7 +83,9 @@ subroutine DAMASK_interface_init() !-------------------------------------------------------------------------------------------------- ! PETSc Init #ifdef _OPENMP - call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr) ! in case of OpenMP, don't rely on PETScInitialize doing MPI init + ! If openMP is enabled, check if the MPI libary supports it and initialize accordingly. + ! Otherwise, the first call to PETSc will do the initialization. + call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr) if (threadLevel0) & - write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) - write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) - write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg) - write(6,'(a,a)') ' Working directory: ', trim(getSolverWorkingDirectoryName()) - write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) - write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) - write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) - if (SpectralRestartInc > 1_pInt) & - write(6,'(a,i6.6)') ' Restart at increment: ', spectralRestartInc - write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile + write(6,'(a,a)') ' Host name: ', trim(hostName) + write(6,'(a,a)') ' User name: ', trim(userName) + write(6,'(a,a)') ' Command line call: ', trim(commandLine) + if (len(trim(workingDirArg)) > 0) & + write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) + write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) + write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg) + write(6,'(a,a)') ' Working directory: ', trim(getSolverWorkingDirectoryName()) + write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) + write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) + write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) + if (SpectralRestartInc > 0_pInt) & + write(6,'(a,i6.6)') ' Restart from increment: ', spectralRestartInc + write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile end subroutine DAMASK_interface_init diff --git a/src/spectral_mech_AL.f90 b/src/spectral_mech_AL.f90 index 4695d4faa..27a1def0c 100644 --- a/src/spectral_mech_AL.f90 +++ b/src/spectral_mech_AL.f90 @@ -43,17 +43,20 @@ module spectral_mech_AL !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. real(pReal), private, dimension(3,3) :: & - F_aimDot, & !< assumed rate of average deformation gradient + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient F_av = 0.0_pReal, & !< average incompatible def grad field P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress P_avLastEval = 0.0_pReal !< average 1st Piola--Kirchhoff stress last call of CPFEM_general - character(len=1024), private :: incInfo !< time and increment information + + character(len=1024), private :: incInfo !< time and increment information + real(pReal), private, dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness S = 0.0_pReal, & !< current compliance (filled up with zeros) C_scale = 0.0_pReal, & S_scale = 0.0_pReal @@ -62,7 +65,7 @@ module spectral_mech_AL err_BC, & !< deviation from stress BC err_curl, & !< RMS of curl of F err_div !< RMS of div of P - logical, private :: ForwardData + integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment @@ -80,11 +83,11 @@ module spectral_mech_AL contains !-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc) !-------------------------------------------------------------------------------------------------- subroutine AL_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -102,12 +105,15 @@ subroutine AL_init use numerics, only: & worldrank, & worldsize + use homogenization, only: & + materialpoint_F0 use DAMASK_interface, only: & getSolverJobName use spectral_utilities, only: & Utilities_constitutiveResponse, & Utilities_updateGamma, & - Utilities_updateIPcoords + Utilities_updateIPcoords, & + wgt use mesh, only: & grid, & grid3 @@ -120,7 +126,11 @@ subroutine AL_init temp33_Real = 0.0_pReal PetscErrorCode :: ierr - PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda + PetscScalar, pointer, dimension(:,:,:,:) :: & + FandF_lambda, & ! overall pointer to solution data + F, & ! specific (sub)pointer + F_lambda ! specific (sub)pointer + integer(pInt), dimension(:), allocatable :: localK integer(pInt) :: proc character(len=1024) :: rankStr @@ -160,83 +170,74 @@ subroutine AL_init DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point grid(1),grid(2),grid(3), & ! global grid 1 , 1, worldsize, & - 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) + 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) grid(1),grid(2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) - call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,PETSC_NULL_OBJECT,ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da + call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESSetConvergenceTest(snes,AL_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" CHKERRQ(ierr) - call SNESSetConvergenceTest(snes,AL_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) - CHKERRQ(ierr) - call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) + call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data - F => xx_psc(0:8,:,:,:) - F_lambda => xx_psc(9:17,:,:,:) - restart: if (restartInc > 1_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & - write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading values of increment ', restartInc - 1_pInt, ' from file' - flush(6) + call DMDAVecGetArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr) ! places pointer on PETSc data + F => FandF_lambda( 0: 8,:,:,:) + F_lambda => FandF_lambda( 9:17,:,:,:) + + restart: if (restartInc > 0_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading values of increment ', restartInc, ' from file' + flush(6) + endif write(rankStr,'(a1,i0)')'_',worldrank call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) + read (777,rec=1) F; close (777) call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F_lastInc - close (777) + read (777,rec=1) F_lastInc; close (777) call IO_read_realFile(777,'F_lambda'//trim(rankStr),trim(getSolverJobName()),size(F_lambda)) - read (777,rec=1) F_lambda - close (777) - call IO_read_realFile(777,'F_lambda_lastInc'//trim(rankStr),& - trim(getSolverJobName()),size(F_lambda_lastInc)) - read (777,rec=1) F_lambda_lastInc - close (777) - call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim)) - read (777,rec=1) F_aim - close (777) - call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc)) - read (777,rec=1) F_aim_lastInc - close (777) - call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) - read (777,rec=1) f_aimDot - close (777) - elseif (restartInc == 1_pInt) then restart + read (777,rec=1) F_lambda; close (777) + call IO_read_realFile(777,'F_lambda_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lambda_lastInc)) + read (777,rec=1) F_lambda_lastInc; close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot)) + read (777,rec=1) F_aimDot; close (777) + F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F + F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc + elseif (restartInc == 0_pInt) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) F_lambda = F F_lambda_lastInc = F_lastInc endif restart + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & - reshape(F,shape(F_lastInc)), 0.0_pReal, math_I3) - + call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + reshape(F,shape(F_lastInc)), & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition nullify(F) nullify(F_lambda) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr) ! write data back to PETSc - restartRead: if (restartInc > 1_pInt) then + restartRead: if (restartInc > 0_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & - write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading more values of increment', restartInc - 1_pInt, 'from file' + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading more values of increment ', restartInc, ' from file' flush(6) call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) - read (777,rec=1) C_volAvg - close (777) + read (777,rec=1) C_volAvg; close (777) call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) - read (777,rec=1) C_volAvgLastInc - close (777) + read (777,rec=1) C_volAvgLastInc; close (777) call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) - read (777,rec=1) C_minMaxAvg - close (777) + read (777,rec=1) C_minMaxAvg; close (777) endif restartRead - call Utilities_updateGamma(C_minMaxAvg,.True.) + call Utilities_updateGamma(C_minMaxAvg,.true.) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) @@ -246,8 +247,7 @@ end subroutine AL_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the AL scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function & - AL_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) +type(tSolutionState) function AL_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) use IO, only: & IO_error use numerics, only: & @@ -266,13 +266,13 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! input data for solution - real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment - type(tBoundaryCondition), intent(in) :: & - stress_BC character(len=*), intent(in) :: & incInfoIn + real(pReal), intent(in) :: & + timeinc, & !< increment time for current solution + timeinc_old !< increment time of last successful increment + type(tBoundaryCondition), intent(in) :: & + stress_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC !-------------------------------------------------------------------------------------------------- @@ -305,18 +305,17 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) - CHKERRQ(ierr) + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(snes,reason,ierr) - CHKERRQ(ierr) + call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) + + AL_solution%converged = reason > 0 + AL_solution%iterationsNeeded = totalIter AL_solution%termIll = terminallyIll terminallyIll = .false. - if (reason == -4) call IO_error(893_pInt) - if (reason < 1) AL_solution%converged = .false. - AL_solution%iterationsNeeded = totalIter + if (reason == -4) call IO_error(893_pInt) ! MPI error end function AL_solution @@ -331,8 +330,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) polarAlpha, & polarBeta use mesh, only: & - grid3, & - grid + grid, & + grid3 use IO, only: & IO_intOut use math, only: & @@ -341,6 +340,10 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) math_mul3333xx33, & math_invSym3333, & math_mul33x33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation use spectral_utilities, only: & wgt, & tensorField_real, & @@ -350,25 +353,17 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) Utilities_constitutiveResponse, & Utilities_divergenceRMS, & Utilities_curlRMS - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation use homogenization, only: & materialpoint_dPdF use FEsolving, only: & terminallyIll implicit none -!-------------------------------------------------------------------------------------------------- -! strange syntax in the next line because otherwise macros expand beyond 132 character limit - DMDALocalInfo, dimension(& - DMDA_LOCAL_INFO_SIZE) :: & - in + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in PetscScalar, & - target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal + target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal !< what is this? PetscScalar, & - target, dimension(3,3,2, X_RANGE, Y_RANGE, Z_RANGE), intent(out) :: f_scal + target, dimension(3,3,2, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: f_scal !< what is this? PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & F_lambda, & @@ -386,33 +381,33 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) SNESGetNumberFunctionEvals, & SNESGetIterationNumber - F => x_scal(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) - F_lambda => x_scal(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) - residual_F => f_scal(1:3,1:3,1,& - X_RANGE,Y_RANGE,Z_RANGE) + F => x_scal(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_lambda => x_scal(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => f_scal(1:3,1:3,1,& + X_RANGE, Y_RANGE, Z_RANGE) residual_F_lambda => f_scal(1:3,1:3,2,& - X_RANGE,Y_RANGE,Z_RANGE) - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + X_RANGE, Y_RANGE, Z_RANGE) F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - newIteration: if(totalIter <= PETScIter) then + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment !-------------------------------------------------------------------------------------------------- -! report begin of new iteration +! begin of new iteration + newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1_pInt - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & - ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & - math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & - math_transpose33(F_aim) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim (lab) =', math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim =', math_transpose33(F_aim) flush(6) endif newIteration @@ -424,7 +419,6 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3)) - enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- @@ -434,8 +428,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) call utilities_FFTtensorBackward() !-------------------------------------------------------------------------------------------------- -! constructing residual - residual_F_lambda = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) +! constructing F_lambda residual + residual_F_lambda = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !< eq (16) in doi: 10.1016/j.ijplas.2014.02.006 !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -444,15 +438,13 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) call Utilities_constitutiveResponse(residual_F,P_av,C_volAvg,C_minMaxAvg, & F - residual_F_lambda/polarBeta,params%timeinc, params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) - ForwardData = .False. !-------------------------------------------------------------------------------------------------- ! calculate divergence tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise call utilities_FFTtensorForward() - err_div = Utilities_divergenceRMS() - call utilities_FFTtensorBackward() + err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -463,7 +455,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) residual_F(1:3,1:3,i,j,k) - & math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3))) & - + residual_F_lambda(1:3,1:3,i,j,k) + + residual_F_lambda(1:3,1:3,i,j,k) !< eq (16) in doi: 10.1016/j.ijplas.2014.02.006 enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- @@ -472,8 +464,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F call utilities_FFTtensorForward() err_curl = Utilities_curlRMS() - call utilities_FFTtensorBackward() - + + nullify(F) + nullify(F_lambda) + nullify(residual_F) + nullify(residual_F_lambda) end subroutine AL_formResidual @@ -488,8 +483,8 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr err_div_tolAbs, & err_curl_tolRel, & err_curl_tolAbs, & - err_stress_tolAbs, & - err_stress_tolRel + err_stress_tolRel, & + err_stress_tolAbs use math, only: & math_mul3333xx33 use FEsolving, only: & @@ -508,24 +503,24 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr real(pReal) :: & curlTol, & divTol, & - BC_tol + BCTol !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & - mask_stress *(P_av - params%stress_BC))) ! mask = 0.0 for no bc + err_BC = maxval(abs((1.0_pReal-mask_stress) * math_mul3333xx33(C_scale,F_aim-F_av) + & + mask_stress * (P_av-params%stress_BC))) ! mask = 0.0 for no bc !-------------------------------------------------------------------------------------------------- ! error calculation - curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_tolAbs) - divTol = max(maxval(abs(P_av)) *err_div_tolRel,err_div_tolAbs) - BC_tol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) converged: if ((totalIter >= itmin .and. & - all([ err_div/divTol, & + all([ err_div /divTol, & err_curl/curlTol, & - err_BC/BC_tol ] < 1.0_pReal)) & + err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then converged @@ -537,12 +532,12 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr !-------------------------------------------------------------------------------------------------- ! report write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' write(6,'(/,a)') ' ===========================================================================' flush(6) @@ -550,132 +545,142 @@ end subroutine AL_converged !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine +!> @details find new boundary conditions and best F estimate for end of current timestep +!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates !-------------------------------------------------------------------------------------------------- subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) - use math, only: & - math_mul33x33, & - math_mul3333xx33, & - math_transpose33, & - math_rotate_backward33 - use numerics, only: & - worldrank - use mesh, only: & - grid3, & - grid - use spectral_utilities, only: & - Utilities_calculateRate, & - Utilities_forwardField, & - Utilities_updateIPcoords, & - tBoundaryCondition, & - cutBack - use IO, only: & - IO_write_JobRealFile - use FEsolving, only: & - restartWrite + use math, only: & + math_mul33x33, & + math_mul3333xx33, & + math_transpose33, & + math_rotate_backward33 + use numerics, only: & + worldrank + use homogenization, only: & + materialpoint_F0 + use mesh, only: & + grid, & + grid3 + use CPFEM2, only: & + CPFEM_age + use spectral_utilities, only: & + Utilities_calculateRate, & + Utilities_forwardField, & + Utilities_updateIPcoords, & + tBoundaryCondition, & + cutBack + use IO, only: & + IO_write_JobRealFile + use FEsolving, only: & + restartWrite - implicit none - real(pReal), intent(in) :: & - timeinc_old, & - timeinc, & - loadCaseTime !< remaining time of current load case - type(tBoundaryCondition), intent(in) :: & - stress_BC, & - deformation_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC - logical, intent(in) :: & - guess - PetscErrorCode :: ierr - PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda - integer(pInt) :: i, j, k - real(pReal), dimension(3,3) :: F_lambda33 - character(len=1024) :: rankStr + implicit none + logical, intent(in) :: & + guess + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + stress_BC, & + deformation_BC + real(pReal), dimension(3,3), intent(in) ::& + rotation_BC + PetscErrorCode :: ierr + PetscScalar, dimension(:,:,:,:), pointer :: FandF_lambda, F, F_lambda + integer(pInt) :: i, j, k + real(pReal), dimension(3,3) :: F_lambda33 + character(len=32) :: rankStr !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) - F => xx_psc(0:8,:,:,:) - F_lambda => xx_psc(9:17,:,:,:) - if (restartWrite) then - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - write(rankStr,'(a1,i0)')'_',worldrank - call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file - write (777,rec=1) F - close (777) - call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file - write (777,rec=1) F_lastInc - close (777) - call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file - write (777,rec=1) F_lambda - close (777) - call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file - write (777,rec=1) F_lambda_lastInc - close (777) - if (worldrank == 0_pInt) then - call IO_write_jobRealFile(777,'F_aim',size(F_aim)) - write (777,rec=1) F_aim - close(777) - call IO_write_jobRealFile(777,'F_aim_lastInc',size(F_aim_lastInc)) - write (777,rec=1) F_aim_lastInc - close(777) - call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) - write (777,rec=1) F_aimDot - close(777) - call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg)) - write (777,rec=1) C_volAvg - close(777) - call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) - write (777,rec=1) C_volAvgLastInc - close(777) - endif - endif + call DMDAVecGetArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr) + F => FandF_lambda( 0: 8,:,:,:) + F_lambda => FandF_lambda( 9:17,:,:,:) - call utilities_updateIPcoords(F) + if (cutBack) then + C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? + C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? + else + !-------------------------------------------------------------------------------------------------- + ! restart information for spectral solver + if (restartWrite) then ! QUESTION: where is this logical properly set? + write(6,'(/,a)') ' writing converged results for restart' + flush(6) - if (cutBack) then - F_aim = F_aim_lastInc - F_lambda = reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid3]) - F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) - C_volAvg = C_volAvgLastInc - else - ForwardData = .True. - C_volAvgLastInc = C_volAvg -!-------------------------------------------------------------------------------------------------- -! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F - f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim) - elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed - f_aimDot = deformation_BC%maskFloat * deformation_BC%values - elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed - f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime - endif - if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old - F_aim_lastInc = F_aim + if (worldrank == 0_pInt) then + call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg)) + write (777,rec=1) C_volAvg; close(777) + call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) + write (777,rec=1) C_volAvgLastInc; close(777) + ! call IO_write_jobRealFile(777,'C_minMaxAvg',size(C_volAvg)) + ! write (777,rec=1) C_minMaxAvg; close(777) + ! call IO_write_jobRealFile(777,'C_minMaxAvgLastInc',size(C_volAvgLastInc)) + ! write (777,rec=1) C_minMaxAvgLastInc; close(777) + call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) + write (777,rec=1) F_aimDot; close(777) + endif + + write(rankStr,'(a1,i0)')'_',worldrank + call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file + write (777,rec=1) F; close (777) + call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_lastInc; close (777) + call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file + write (777,rec=1) F_lambda; close (777) + call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_lambda_lastInc; close (777) + endif + + call CPFEM_age() ! age state and kinematics + call utilities_updateIPcoords(F) + + C_volAvgLastInc = C_volAvg + C_minMaxAvgLastInc = C_minMaxAvg + + if (guess) then ! QUESTION: better with a = L ? x:y + F_aimDot = stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old ! initialize with correction based on last inc + else + F_aimDot = 0.0_pReal + endif + F_aim_lastInc = F_aim + !-------------------------------------------------------------------------------------------------- + ! calculate rate for aim + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc) + elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * deformation_BC%values + elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + endif + + + Fdot = Utilities_calculateRate(guess, & + F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & + math_rotate_backward33(F_aimDot,rotation_BC)) + F_lambdaDot = Utilities_calculateRate(guess, & + F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid3]), timeinc_old, & + math_rotate_backward33(F_aimDot,rotation_BC)) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward + F_lambda_lastInc = reshape(F_lambda, [3,3,grid(1),grid(2),grid3]) ! winding F_lambda forward + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + endif !-------------------------------------------------------------------------------------------------- -! update coordinates and rate and forward last inc - call utilities_updateIPcoords(F) - Fdot = Utilities_calculateRate(guess, & - F_lastInc, reshape(F, [3,3,grid(1),grid(2),grid3]), timeinc_old, & - math_rotate_backward33(f_aimDot,rotation_BC)) - F_lambdaDot = Utilities_calculateRate(guess, & - F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid3]), timeinc_old, & - math_rotate_backward33(f_aimDot,rotation_BC)) - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) - F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid3]) - endif +! update average and local deformation gradients + F_aim = F_aim_lastInc + F_aimDot * timeinc - F_aim = F_aim + f_aimDot * timeinc - -!-------------------------------------------------------------------------------------------------- -! update local deformation gradient - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim - math_rotate_backward33(F_aim,rotation_BC)), & - [9,grid(1),grid(2),grid3]) - F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & - [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition - if (.not. guess) then ! large strain forwarding - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + math_rotate_backward33(F_aim,rotation_BC)),& + [9,grid(1),grid(2),grid3]) + if (guess) then + F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & + [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + else + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) F_lambda33 = reshape(F_lambda(1:9,i,j,k),[3,3]) F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & math_mul3333xx33(C_scale,& @@ -683,9 +688,12 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stre F_lambda33) -math_I3))*0.5_pReal)& + math_I3 F_lambda(1:9,i,j,k) = reshape(F_lambda33,[9]) - enddo; enddo; enddo - endif - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + enddo; enddo; enddo + endif + + nullify(F) + nullify(F_lambda) + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr) end subroutine AL_forward diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index 0b3409b52..171eeacad 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -38,21 +38,27 @@ module spectral_mech_basic !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. real(pReal), private, dimension(3,3) :: & - F_aim = math_I3, & - F_aim_lastInc = math_I3, & - P_av = 0.0_pReal, & - F_aimDot = 0.0_pReal - character(len=1024), private :: incInfo + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient + F_aim = math_I3, & !< current prescribed deformation gradient + F_aim_lastInc = math_I3, & !< previous average deformation gradient + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress + + character(len=1024), private :: incInfo !< time and increment information + real(pReal), private, dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness S = 0.0_pReal !< current compliance (filled up with zeros) - real(pReal), private :: err_stress, err_div - logical, private :: ForwardData + + real(pReal), private :: & + err_BC, & !< deviation from stress BC + err_div !< RMS of div of P + integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment + real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal public :: & @@ -72,7 +78,7 @@ contains !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine basicPETSc_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -134,8 +140,8 @@ subroutine basicPETSc_init !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -154,11 +160,10 @@ subroutine basicPETSc_init grid(1),grid(2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da call SNESSetConvergenceTest(snes,BasicPETSC_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments @@ -167,10 +172,10 @@ subroutine basicPETSc_init ! init fields call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with - restart: if (restartInc > 1_pInt) then + restart: if (restartInc > 0_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then - write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading values of increment ', restartInc - 1_pInt, ' from file' + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading values of increment ', restartInc, ' from file' flush(6) endif write(rankStr,'(a1,i0)')'_',worldrank @@ -182,24 +187,24 @@ subroutine basicPETSc_init read (777,rec=1) F_aimDot; close (777) F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc - elseif (restartInc == 1_pInt) then restart + elseif (restartInc == 0_pInt) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) endif restart materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P, temp33_Real, C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F 0.0_pReal, & ! time increment math_I3) ! no rotation of boundary condition call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc ! QUESTION: why not writing back right after reading (l.189)? - restartRead: if (restartInc > 1_pInt) then ! QUESTION: are those values not calc'ed by constitutiveResponse? why reading from file? - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & - write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading more values of increment', restartInc-1_pInt, 'from file' + restartRead: if (restartInc > 0_pInt) then ! QUESTION: are those values not calc'ed by constitutiveResponse? why reading from file? + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0 .and. worldrank == 0_pInt) & + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading more values of increment ', restartInc, ' from file' flush(6) call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) read (777,rec=1) C_volAvg; close (777) @@ -208,11 +213,11 @@ subroutine basicPETSc_init call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) read (777,rec=1) C_minMaxAvg; close (777) endif restartRead - - call Utilities_updateGamma(C_minmaxAvg,.true.) + + call Utilities_updateGamma(C_minMaxAvg,.true.) end subroutine basicPETSc_init - + !-------------------------------------------------------------------------------------------------- !> @brief solution for the Basic PETSC scheme with internal iterations !-------------------------------------------------------------------------------------------------- @@ -256,7 +261,7 @@ type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old, !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite) + if (update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite) !-------------------------------------------------------------------------------------------------- @@ -269,8 +274,7 @@ type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old, !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) - CHKERRQ(ierr) + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence @@ -282,7 +286,6 @@ type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old, terminallyIll = .false. if (reason == -4) call IO_error(893_pInt) ! MPI error - end function BasicPETSc_solution @@ -298,7 +301,6 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) grid3 use math, only: & math_rotate_backward33, & - math_transpose33, & math_mul3333xx33 use debug, only: & debug_level, & @@ -346,15 +348,15 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', math_transpose33(F_aim) + ' deformation gradient aim =', transpose(F_aim) flush(6) endif newIteration !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(f_scal,P_av,C_volAvg,C_minmaxAvg, & + call Utilities_constitutiveResponse(f_scal,P_av,C_volAvg,C_minMaxAvg, & x_scal,params%timeinc, params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) @@ -362,7 +364,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) ! stress BC handling deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) F_aim = F_aim - deltaF_aim - err_stress = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + err_BC = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme @@ -406,14 +408,14 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du PetscErrorCode :: ierr real(pReal) :: & divTol, & - stressTol + BCTol - divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs) - stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) + divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) converged: if ((totalIter >= itmin .and. & all([ err_div/divTol, & - err_stress/stressTol ] < 1.0_pReal)) & + err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then converged @@ -426,9 +428,9 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du ! report write(6,'(1/,a)') ' ... reporting .............................................................' write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' - write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & - err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' + err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' write(6,'(/,a)') ' ===========================================================================' flush(6) @@ -445,8 +447,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation math_rotate_backward33 use numerics, only: & worldrank - use homogenization, only: & - materialpoint_F0 + use homogenization, only: & + materialpoint_F0 use mesh, only: & grid, & grid3 @@ -476,7 +478,7 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation real(pReal), dimension(3,3), intent(in) ::& rotation_BC PetscErrorCode :: ierr - PetscScalar, pointer :: F(:,:,:,:) + PetscScalar, dimension(:,:,:,:), pointer :: F character(len=32) :: rankStr @@ -497,10 +499,6 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation write (777,rec=1) C_volAvg; close(777) call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) write (777,rec=1) C_volAvgLastInc; close(777) - call IO_write_jobRealFile(777,'C_minMaxAvg',size(C_volAvg)) - write (777,rec=1) C_minMaxAvg; close(777) - call IO_write_jobRealFile(777,'C_minMaxAvgLastInc',size(C_volAvgLastInc)) - write (777,rec=1) C_minMaxAvgLastInc; close(777) call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) write (777,rec=1) F_aimDot; close(777) endif @@ -517,13 +515,10 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - - if (guess) then ! QUESTION: better with a = L ? x:y - F_aimDot = stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old ! initialize with correction based on last inc - else - F_aimDot = 0.0_pReal - endif + + F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) F_aim_lastInc = F_aim + !-------------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F @@ -551,7 +546,7 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) - + end subroutine BasicPETSc_forward !-------------------------------------------------------------------------------------------------- diff --git a/src/spectral_mech_Polarisation.f90 b/src/spectral_mech_Polarisation.f90 index fc65f14cf..acd713c70 100644 --- a/src/spectral_mech_Polarisation.f90 +++ b/src/spectral_mech_Polarisation.f90 @@ -43,7 +43,7 @@ module spectral_mech_Polarisation !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. real(pReal), private, dimension(3,3) :: & - F_aimDot, & !< assumed rate of average deformation gradient + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient F_av = 0.0_pReal, & !< average incompatible def grad field @@ -54,6 +54,7 @@ module spectral_mech_Polarisation C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness S = 0.0_pReal, & !< current compliance (filled up with zeros) C_scale = 0.0_pReal, & S_scale = 0.0_pReal @@ -62,7 +63,7 @@ module spectral_mech_Polarisation err_BC, & !< deviation from stress BC err_curl, & !< RMS of curl of F err_div !< RMS of div of P - logical, private :: ForwardData + integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment @@ -80,11 +81,11 @@ module spectral_mech_Polarisation contains !-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc) !-------------------------------------------------------------------------------------------------- subroutine Polarisation_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -102,12 +103,15 @@ subroutine Polarisation_init use numerics, only: & worldrank, & worldsize + use homogenization, only: & + materialpoint_F0 use DAMASK_interface, only: & getSolverJobName use spectral_utilities, only: & Utilities_constitutiveResponse, & Utilities_updateGamma, & - Utilities_updateIPcoords + Utilities_updateIPcoords, & + wgt use mesh, only: & grid, & grid3 @@ -120,7 +124,11 @@ subroutine Polarisation_init temp33_Real = 0.0_pReal PetscErrorCode :: ierr - PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_tau + PetscScalar, pointer, dimension(:,:,:,:) :: & + FandF_tau, & ! overall pointer to solution data + F, & ! specific (sub)pointer + F_tau ! specific (sub)pointer + integer(pInt), dimension(:), allocatable :: localK integer(pInt) :: proc character(len=1024) :: rankStr @@ -164,78 +172,69 @@ subroutine Polarisation_init grid(1),grid(2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) - call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,PETSC_NULL_OBJECT,ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da + call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESSetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" CHKERRQ(ierr) - call SNESSetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) - CHKERRQ(ierr) - call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) + call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data - F => xx_psc(0:8,:,:,:) - F_tau => xx_psc(9:17,:,:,:) - restart: if (restartInc > 1_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & - write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading values of increment ', restartInc - 1_pInt, ' from file' - flush(6) + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data + F => FandF_tau( 0: 8,:,:,:) + F_tau => FandF_tau( 9:17,:,:,:) + restart: if (restartInc > 0_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading values of increment ', restartInc, ' from file' + flush(6) + endif write(rankStr,'(a1,i0)')'_',worldrank call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) + read (777,rec=1) F; close (777) call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F_lastInc - close (777) + read (777,rec=1) F_lastInc; close (777) call IO_read_realFile(777,'F_tau'//trim(rankStr),trim(getSolverJobName()),size(F_tau)) - read (777,rec=1) F_tau - close (777) - call IO_read_realFile(777,'F_tau_lastInc'//trim(rankStr),& - trim(getSolverJobName()),size(F_tau_lastInc)) - read (777,rec=1) F_tau_lastInc - close (777) - call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim)) - read (777,rec=1) F_aim - close (777) - call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc)) - read (777,rec=1) F_aim_lastInc - close (777) - call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) - read (777,rec=1) f_aimDot - close (777) - elseif (restartInc == 1_pInt) then restart - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + read (777,rec=1) F_tau; close (777) + call IO_read_realFile(777,'F_tau_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_tau_lastInc)) + read (777,rec=1) F_tau_lastInc; close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot)) + read (777,rec=1) F_aimDot; close (777) + F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F + F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc + elseif (restartInc == 0_pInt) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) - F_tau = 2.0_pReal* F + F_tau = 2.0_pReal*F F_tau_lastInc = 2.0_pReal*F_lastInc endif restart + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & - reshape(F,shape(F_lastInc)),0.0_pReal,math_I3) + call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + reshape(F,shape(F_lastInc)), & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition nullify(F) nullify(F_tau) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! write data back to PETSc - restartRead: if (restartInc > 1_pInt) then + restartRead: if (restartInc > 0_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & - write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading more values of increment', restartInc - 1_pInt, 'from file' + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading more values of increment ', restartInc, ' from file' flush(6) call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) - read (777,rec=1) C_volAvg - close (777) + read (777,rec=1) C_volAvg; close (777) call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) - read (777,rec=1) C_volAvgLastInc - close (777) + read (777,rec=1) C_volAvgLastInc; close (777) call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) - read (777,rec=1) C_minMaxAvg - close (777) + read (777,rec=1) C_minMaxAvg; close (777) endif restartRead - call Utilities_updateGamma(C_minMaxAvg,.True.) + call Utilities_updateGamma(C_minMaxAvg,.true.) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) @@ -245,8 +244,7 @@ end subroutine Polarisation_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function & - Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) +type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) use IO, only: & IO_error use numerics, only: & @@ -265,13 +263,13 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! input data for solution - real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment - type(tBoundaryCondition), intent(in) :: & - stress_BC character(len=*), intent(in) :: & incInfoIn + real(pReal), intent(in) :: & + timeinc, & !< increment time for current solution + timeinc_old !< increment time of last successful increment + type(tBoundaryCondition), intent(in) :: & + stress_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC !-------------------------------------------------------------------------------------------------- @@ -304,18 +302,17 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) - CHKERRQ(ierr) + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(snes,reason,ierr) - CHKERRQ(ierr) + call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) + + Polarisation_solution%converged = reason > 0 + Polarisation_solution%iterationsNeeded = totalIter Polarisation_solution%termIll = terminallyIll terminallyIll = .false. - if (reason == -4) call IO_error(893_pInt) - if (reason < 1) Polarisation_solution%converged = .false. - Polarisation_solution%iterationsNeeded = totalIter + if (reason == -4) call IO_error(893_pInt) ! MPI error end function Polarisation_solution @@ -330,16 +327,19 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) polarAlpha, & polarBeta use mesh, only: & - grid3, & - grid + grid, & + grid3 use IO, only: & IO_intOut use math, only: & math_rotate_backward33, & - math_transpose33, & math_mul3333xx33, & math_invSym3333, & math_mul33x33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation use spectral_utilities, only: & wgt, & tensorField_real, & @@ -349,21 +349,13 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) Utilities_constitutiveResponse, & Utilities_divergenceRMS, & Utilities_curlRMS - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation use homogenization, only: & materialpoint_dPdF use FEsolving, only: & terminallyIll implicit none -!-------------------------------------------------------------------------------------------------- -! strange syntax in the next line because otherwise macros expand beyond 132 character limit - DMDALocalInfo, dimension(& - DMDA_LOCAL_INFO_SIZE) :: & - in + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in PetscScalar, & target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal PetscScalar, & @@ -385,33 +377,33 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) SNESGetNumberFunctionEvals, & SNESGetIterationNumber - F => x_scal(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) - F_tau => x_scal(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) - residual_F => f_scal(1:3,1:3,1,& - X_RANGE,Y_RANGE,Z_RANGE) - residual_F_tau => f_scal(1:3,1:3,2,& - X_RANGE,Y_RANGE,Z_RANGE) - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + F => x_scal(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_tau => x_scal(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => f_scal(1:3,1:3,1,& + X_RANGE, Y_RANGE, Z_RANGE) + residual_F_tau => f_scal(1:3,1:3,2,& + X_RANGE, Y_RANGE, Z_RANGE) F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - newIteration: if(totalIter <= PETScIter) then + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment !-------------------------------------------------------------------------------------------------- -! report begin of new iteration +! begin of new iteration + newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1_pInt - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & - ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & - math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & - math_transpose33(F_aim) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim =', transpose(F_aim) flush(6) endif newIteration @@ -441,15 +433,13 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) call Utilities_constitutiveResponse(residual_F,P_av,C_volAvg,C_minMaxAvg, & F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) - ForwardData = .False. !-------------------------------------------------------------------------------------------------- ! calculate divergence tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise call utilities_FFTtensorForward() - err_div = Utilities_divergenceRMS() - call utilities_FFTtensorBackward() + err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -469,8 +459,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F call utilities_FFTtensorForward() err_curl = Utilities_curlRMS() - call utilities_FFTtensorBackward() - + + nullify(F) + nullify(F_tau) + nullify(residual_F) + nullify(residual_F_tau) end subroutine Polarisation_formResidual @@ -485,8 +478,8 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, err_div_tolAbs, & err_curl_tolRel, & err_curl_tolAbs, & - err_stress_tolAbs, & - err_stress_tolRel + err_stress_tolRel, & + err_stress_tolAbs use math, only: & math_mul3333xx33 use FEsolving, only: & @@ -505,24 +498,24 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, real(pReal) :: & curlTol, & divTol, & - BC_tol + BCTol !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & - mask_stress *(P_av - params%stress_BC))) ! mask = 0.0 for no bc + err_BC = maxval(abs((1.0_pReal-mask_stress) * math_mul3333xx33(C_scale,F_aim-F_av) + & + mask_stress * (P_av-params%stress_BC))) ! mask = 0.0 for no bc !-------------------------------------------------------------------------------------------------- ! error calculation - curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_tolAbs) - divTol = max(maxval(abs(P_av)) *err_div_tolRel,err_div_tolAbs) - BC_tol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) converged: if ((totalIter >= itmin .and. & - all([ err_div/divTol, & + all([ err_div /divTol, & err_curl/curlTol, & - err_BC/BC_tol ] < 1.0_pReal)) & + err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then converged @@ -534,12 +527,12 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, !-------------------------------------------------------------------------------------------------- ! report write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' write(6,'(/,a)') ' ===========================================================================' flush(6) @@ -547,142 +540,146 @@ end subroutine Polarisation_converged !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine +!> @details find new boundary conditions and best F estimate for end of current timestep +!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates !-------------------------------------------------------------------------------------------------- subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) use math, only: & math_mul33x33, & math_mul3333xx33, & - math_transpose33, & math_rotate_backward33 use numerics, only: & worldrank + use homogenization, only: & + materialpoint_F0 use mesh, only: & - grid3, & - grid - use spectral_utilities, only: & - Utilities_calculateRate, & - Utilities_forwardField, & - Utilities_updateIPcoords, & - tBoundaryCondition, & - cutBack - use IO, only: & - IO_write_JobRealFile - use FEsolving, only: & - restartWrite + grid, & + grid3 + use CPFEM2, only: & + CPFEM_age + use spectral_utilities, only: & + Utilities_calculateRate, & + Utilities_forwardField, & + Utilities_updateIPcoords, & + tBoundaryCondition, & + cutBack + use IO, only: & + IO_write_JobRealFile + use FEsolving, only: & + restartWrite implicit none - real(pReal), intent(in) :: & - timeinc_old, & - timeinc, & - loadCaseTime !< remaining time of current load case - type(tBoundaryCondition), intent(in) :: & - stress_BC, & - deformation_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC - logical, intent(in) :: & - guess - PetscErrorCode :: ierr - PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_tau - integer(pInt) :: i, j, k - real(pReal), dimension(3,3) :: F_lambda33 - character(len=1024) :: rankStr + logical, intent(in) :: & + guess + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + stress_BC, & + deformation_BC + real(pReal), dimension(3,3), intent(in) ::& + rotation_BC + PetscErrorCode :: ierr + PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau + integer(pInt) :: i, j, k + real(pReal), dimension(3,3) :: F_lambda33 + character(len=32) :: rankStr !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) - F => xx_psc(0:8,:,:,:) - F_tau => xx_psc(9:17,:,:,:) - if (restartWrite) then - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - write(rankStr,'(a1,i0)')'_',worldrank - call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file - write (777,rec=1) F - close (777) - call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file - write (777,rec=1) F_lastInc - close (777) - call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file - write (777,rec=1) F_tau - close (777) - call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_lastInc field to file - write (777,rec=1) F_tau_lastInc - close (777) - if (worldrank == 0_pInt) then - call IO_write_jobRealFile(777,'F_aim',size(F_aim)) - write (777,rec=1) F_aim - close(777) - call IO_write_jobRealFile(777,'F_aim_lastInc',size(F_aim_lastInc)) - write (777,rec=1) F_aim_lastInc - close (777) - call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) - write (777,rec=1) F_aimDot - close(777) - call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg)) - write (777,rec=1) C_volAvg - close(777) - call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) - write (777,rec=1) C_volAvgLastInc - close(777) - endif - endif + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) + F => FandF_tau( 0: 8,:,:,:) + F_tau => FandF_tau( 9:17,:,:,:) - call utilities_updateIPcoords(F) + if (cutBack) then + C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? + C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? + else + !-------------------------------------------------------------------------------------------------- + ! restart information for spectral solver + if (restartWrite) then ! QUESTION: where is this logical properly set? + write(6,'(/,a)') ' writing converged results for restart' + flush(6) - if (cutBack) then - F_aim = F_aim_lastInc - F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid3]) - F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) - C_volAvg = C_volAvgLastInc - else - ForwardData = .True. - C_volAvgLastInc = C_volAvg -!-------------------------------------------------------------------------------------------------- -! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F - f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim) - elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed - f_aimDot = deformation_BC%maskFloat * deformation_BC%values - elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed - f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime - endif - if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old - F_aim_lastInc = F_aim + if (worldrank == 0_pInt) then + call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg)) + write (777,rec=1) C_volAvg; close(777) + call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) + write (777,rec=1) C_volAvgLastInc; close(777) + call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) + write (777,rec=1) F_aimDot; close(777) + endif + + write(rankStr,'(a1,i0)')'_',worldrank + call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file + write (777,rec=1) F; close (777) + call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_lastInc; close (777) + call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file + write (777,rec=1) F_tau; close (777) + call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_tau_lastInc field to file + write (777,rec=1) F_tau_lastInc; close (777) + endif + + call CPFEM_age() ! age state and kinematics + call utilities_updateIPcoords(F) + + C_volAvgLastInc = C_volAvg + C_minMaxAvgLastInc = C_minMaxAvg + + F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aim_lastInc = F_aim + !-------------------------------------------------------------------------------------------------- + ! calculate rate for aim + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc) + elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * deformation_BC%values + elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + endif + + + Fdot = Utilities_calculateRate(guess, & + F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & + math_rotate_backward33(F_aimDot,rotation_BC)) + F_tauDot = Utilities_calculateRate(guess, & + F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & + math_rotate_backward33(F_aimDot,rotation_BC)) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward + F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau forward + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + endif !-------------------------------------------------------------------------------------------------- -! update coordinates and rate and forward last inc - call utilities_updateIPcoords(F) - Fdot = Utilities_calculateRate(guess, & - F_lastInc, reshape(F, [3,3,grid(1),grid(2),grid3]), timeinc_old, & - math_rotate_backward33( f_aimDot,rotation_BC)) - F_tauDot = Utilities_calculateRate(guess, & - F_tau_lastInc, reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & - math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC)) - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) - F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) - endif +! update average and local deformation gradients + F_aim = F_aim_lastInc + F_aimDot * timeinc - F_aim = F_aim + f_aimDot * timeinc - -!-------------------------------------------------------------------------------------------------- -! update local deformation gradient - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim - math_rotate_backward33(F_aim,rotation_BC)), & - [9,grid(1),grid(2),grid3]) - F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & ! does not have any average value as boundary condition - [9,grid(1),grid(2),grid3]) - if (.not. guess) then ! large strain forwarding + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + math_rotate_backward33(F_aim,rotation_BC)),& + [9,grid(1),grid(2),grid3]) + if (guess) then + F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & + [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + else do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & math_mul3333xx33(C_scale,& - math_mul33x33(math_transpose33(F_lambda33),& - F_lambda33) -math_I3))*0.5_pReal)& + math_mul33x33(transpose(F_lambda33),& + F_lambda33)-math_I3))*0.5_pReal)& + math_I3 F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) enddo; enddo; enddo endif - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + + nullify(F) + nullify(F_tau) + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) end subroutine Polarisation_forward diff --git a/src/spectral_thermal.f90 b/src/spectral_thermal.f90 index 322f12031..f89184543 100644 --- a/src/spectral_thermal.f90 +++ b/src/spectral_thermal.f90 @@ -60,7 +60,7 @@ contains !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine spectral_thermal_init -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index a66fa558e..e3383f3d1 100644 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -70,7 +70,6 @@ module spectral_utilities ! derived types type, public :: tSolutionState !< return type of solution from spectral solver variants logical :: converged = .true. - logical :: regrid = .false. logical :: stagConverged = .true. logical :: termIll = .false. integer(pInt) :: iterationsNeeded = 0_pInt @@ -145,8 +144,7 @@ module spectral_utilities FIELD_UNDEFINED_ID, & FIELD_MECH_ID, & FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID, & - utilities_calcPlasticity + FIELD_DAMAGE_ID private :: & utilities_getFreqDerivative @@ -154,14 +152,14 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, sets debug flags, create plans for FFTW -!> @details Sets the debug levels for general, divergence, restart and FFTW from the biwise coding +!> @details Sets the debug levels for general, divergence, restart, and FFTW from the bitwise coding !> provided by the debug module to logicals. -!> Allocates all fields used by FFTW and create the corresponding plans depending on the debug +!> Allocate all fields used by FFTW and create the corresponding plans depending on the debug !> level chosen. !> Initializes FFTW. !-------------------------------------------------------------------------------------------------- subroutine utilities_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options @@ -378,10 +376,10 @@ end subroutine utilities_init !-------------------------------------------------------------------------------------------------- -!> @brief updates references stiffness and potentially precalculated gamma operator +!> @brief updates reference stiffness and potentially precalculated gamma operator !> @details Sets the current reference stiffness to the stiffness given as an argument. !> If the gamma operator is precalculated, it is calculated with this stiffness. -!> In case of a on-the-fly calculation, only the reference stiffness is updated. +!> In case of an on-the-fly calculation, only the reference stiffness is updated. !> Also writes out the current reference stiffness for restart. !-------------------------------------------------------------------------------------------------- subroutine utilities_updateGamma(C,saveReference) @@ -414,8 +412,7 @@ subroutine utilities_updateGamma(C,saveReference) write(6,'(/,a)') ' writing reference stiffness to file' flush(6) call IO_write_jobRealFile(777,'C_ref',size(C_ref)) - write (777,rec=1) C_ref - close(777) + write (777,rec=1) C_ref; close(777) endif endif @@ -784,7 +781,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) if(debugGeneral) then write(6,'(/,a)') ' ... updating masked compliance ............................................' write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& - transpose(temp99_Real)/1.e9_pReal + transpose(temp99_Real)*1.0e-9_pReal flush(6) endif k = 0_pInt ! calculate reduced stiffness @@ -823,7 +820,6 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) .or. (m/=n .and. abs(sTimesC(m,n)) > 1.0e-12_pReal) ! off-diagonal elements of S*C should be 0 enddo enddo - if (debugGeneral .or. errmatinv) then write(formatString, '(i2)') size_reduced formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' @@ -838,13 +834,11 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) else temp99_real = 0.0_pReal endif - if(debugGeneral) then write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') & - ' Masked Compliance (load) / GPa =', transpose(temp99_Real*1.e-9_pReal) + ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal flush(6) endif - utilities_maskedCompliance = math_Plain99to3333(temp99_Real) end function utilities_maskedCompliance @@ -940,7 +934,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& debug_reset, & debug_info use math, only: & - math_transpose33, & math_rotate_forward33, & math_det33 use mesh, only: & @@ -999,10 +992,10 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) if (debugRotation) & write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& - math_transpose33(P_av)*1.e-6_pReal + transpose(P_av)*1.e-6_pReal P_av = math_rotate_forward33(P_av,rotation_BC) write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& - math_transpose33(P_av)*1.e-6_pReal + transpose(P_av)*1.e-6_pReal flush(6) max_dPdF = 0.0_pReal @@ -1034,125 +1027,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& end subroutine utilities_constitutiveResponse -!-------------------------------------------------------------------------------------------------- -!> @brief calculates yield stress, plastic strain, total strain and their equivalent values -!-------------------------------------------------------------------------------------------------- -subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTotalStrain, & - eqPlasticStrain, plasticWork, rotation_BC) - use crystallite, only: & - crystallite_Fe, & - crystallite_P, & - crystallite_subF - use material, only: & - homogenization_maxNgrains - use mesh, only: & - mesh_maxNips,& - mesh_NcpElems - use math, only: & - math_det33, & - math_inv33, & - math_mul33x33, & - math_trace33, & - math_transpose33, & - math_equivStrain33, & - math_equivStress33, & - math_rotate_forward33, & - math_identity2nd, & - math_crossproduct, & - math_eigenvectorBasisSym, & - math_eigenvectorBasisSym33, & - math_eigenvectorBasisSym33_log, & - math_eigenValuesVectorsSym33 - - implicit none - - real(pReal), intent(inout) :: eqStress, eqPlasticStrain, plasticWork - real(pReal), intent(out) :: eqTotalStrain - real(pReal), dimension(3,3),intent(out) :: yieldStress, plasticStrain - real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame - real(pReal), dimension(3,3) :: cauchy, P_av, F_av, Ve_av !< average - real(pReal), dimension(3) :: Values, S - real(pReal), dimension(3,3) :: Vectors, diag - real(pReal), dimension(3,3) :: & - Vp, F_temp, U, VT, R, V, V_total - real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - Be, Ve, Fe - real(pReal), dimension(15) :: WORK !< previous deformation gradient - integer(pInt) :: INFO, i, j, k, l, ierr - real(pReal) :: wgtm - real(pReal) :: eqStressOld, eqPlasticStrainOld, plasticWorkOld - - external :: dgesvd - - eqStressOld = eqStress - eqPlasticStrainOld = eqPlasticStrain - plasticWorkOld = plasticWork - wgtm = 1.0_pReal/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal) - diag = 0.0_pReal - - P_av = sum(sum(sum(crystallite_P,dim=5),dim=4),dim=3) * wgtm - call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - P_av = math_rotate_forward33(P_av,rotation_BC) - - F_av = sum(sum(sum(crystallite_subF,dim=5),dim=4),dim=3) * wgtm - call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - F_av = math_rotate_forward33(F_av,rotation_BC) - - cauchy = 1.0_pReal/math_det33(F_av)*math_mul33x33(P_av,transpose(F_av)) - yieldStress = cauchy - eqStress = math_equivStress33(cauchy) - - F_temp = F_av - call dgesvd ('A', 'A', 3, 3, F_temp, 3, S, U, 3, VT, 3, WORK, 15, INFO) ! singular value decomposition - - R = math_mul33x33(U, VT) ! rotation of polar decomposition - V = math_mul33x33(F_av,math_inv33(R)) - - call math_eigenValuesVectorsSym33(V,Values,Vectors) - do l = 1_pInt, 3_pInt - if (Values(l) < 0.0_pReal) then - Values(l) = -Values(l) - Vectors(1:3, l) = -Vectors(1:3, l) - endif - Values(l) = log(Values(l)) - diag(l,l) = Values(l) - enddo - if (dot_product(Vectors(1:3,1),Vectors(1:3,2)) /= 0) then - Vectors(1:3,2) = math_crossproduct(Vectors(1:3,3), Vectors(1:3,1)) - Vectors(1:3,2) = Vectors(1:3,2)/sqrt(dot_product(Vectors(1:3,2),Vectors(1:3,2))) - endif - if (dot_product(Vectors(1:3,2),Vectors(1:3,3)) /= 0) then - Vectors(1:3,3) = math_crossproduct(Vectors(1:3,1), Vectors(1:3,2)) - Vectors(1:3,3) = Vectors(1:3,3)/sqrt(dot_product(Vectors(1:3,3),Vectors(1:3,3))) - endif - if (dot_product(Vectors(1:3,3),Vectors(1:3,1)) /= 0) then - Vectors(1:3,1) = math_crossproduct(Vectors(1:3,2), Vectors(1:3,3)) - Vectors(1:3,1) = Vectors(1:3,1)/sqrt(dot_product(Vectors(1:3,1),Vectors(1:3,1))) - endif - - V_total = REAL(math_mul33x33(Vectors, math_mul33x33(diag, transpose(Vectors)))) - eqTotalStrain = math_equivStrain33(V_total) - - do k = 1_pInt, mesh_NcpElems; do j = 1_pInt, mesh_maxNips; do i = 1_pInt,homogenization_maxNgrains - Fe(1:3,1:3,i,j,k) = crystallite_Fe(1:3,1:3,i,j,k) - Fe(1:3,1:3,i,j,k) = math_rotate_forward33(Fe(1:3,1:3,i,j,k),rotation_BC) - Be(1:3,1:3,i,j,k) = math_mul33x33(Fe(1:3,1:3,i,j,k),math_transpose33(Fe(1:3,1:3,i,j,k))) ! elastic part of left Cauchy–Green deformation tensor - Ve(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Be(1:3,1:3,i,j,k)) - enddo; enddo; enddo - - Ve_av = sum(sum(sum(Ve,dim=5),dim=4),dim=3) * wgtm - call MPI_Allreduce(MPI_IN_PLACE,Ve_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - - Vp = V_total - Ve_av - - eqPlasticStrain = math_equivStrain33(Vp) - - plasticStrain = Vp - - plasticWork = plasticWorkOld + 0.5*(eqStressOld + eqStress) * (eqPlasticStrain - eqPlasticStrainOld) - -end subroutine utilities_calcPlasticity - !-------------------------------------------------------------------------------------------------- !> @brief calculates forward rate, either guessing or just add delta/timeinc diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 7f23a81b5..62ffabf9c 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -46,7 +46,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index c55d1d3eb..151eb7aa3 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -47,7 +47,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine thermal_conduction_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index 87e846f12..68e09de8c 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -16,7 +16,7 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- subroutine thermal_isothermal_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/vacancyflux_cahnhilliard.f90 b/src/vacancyflux_cahnhilliard.f90 index 9f6ecd8b0..e40772d11 100644 --- a/src/vacancyflux_cahnhilliard.f90 +++ b/src/vacancyflux_cahnhilliard.f90 @@ -61,7 +61,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine vacancyflux_cahnhilliard_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/vacancyflux_isochempot.f90 b/src/vacancyflux_isochempot.f90 index 8c256467f..f98379eba 100644 --- a/src/vacancyflux_isochempot.f90 +++ b/src/vacancyflux_isochempot.f90 @@ -44,7 +44,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine vacancyflux_isochempot_init(fileUnit) -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options diff --git a/src/vacancyflux_isoconc.f90 b/src/vacancyflux_isoconc.f90 index ad7842e3f..470560206 100644 --- a/src/vacancyflux_isoconc.f90 +++ b/src/vacancyflux_isoconc.f90 @@ -16,7 +16,7 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- subroutine vacancyflux_isoconc_init() -#ifdef __GFORTRAN__ +#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options