Merge remote-tracking branch 'origin/development' into physics-based-hex-interactions

.github/workflows/Fortran.yml

@@ -2,11 +2,12 @@ name: Grid and Mesh Solver
 on: [push]
 
 env:
-  HOMEBREW_NO_ANALYTICS: "ON"     # Make Homebrew installation a little quicker
+  PETSC_VERSION: '3.16.2'
-  HOMEBREW_NO_AUTO_UPDATE: "ON"
+  HOMEBREW_NO_ANALYTICS: 'ON'     # Make Homebrew installation a little quicker
-  HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: "ON"
+  HOMEBREW_NO_AUTO_UPDATE: 'ON'
-  HOMEBREW_NO_GITHUB_API: "ON"
+  HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: 'ON'
-  HOMEBREW_NO_INSTALL_CLEANUP: "ON"
+  HOMEBREW_NO_GITHUB_API: 'ON'
+  HOMEBREW_NO_INSTALL_CLEANUP: 'ON'
 
 jobs:
 
@@ -48,17 +49,17 @@ jobs:
         uses: actions/cache@v2
         with:
           path: download
-          key: petsc-3.16.0.tar.gz
+          key: petsc-${{ env.PETSC_VERSION }}.tar.gz
 
       - name: PETSc - Download
         if: steps.petsc-download.outputs.cache-hit != 'true'
         run: |
-          wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.16.0.tar.gz -P download
+          wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-${PETSC_VERSION}.tar.gz -P download
 
       - name: PETSc - Prepare
         run: |
-          tar -xf download/petsc-3.16.0.tar.gz -C .
+          tar -xf download/petsc-${PETSC_VERSION}.tar.gz -C .
-          export PETSC_DIR=${PWD}/petsc-3.16.0
+          export PETSC_DIR=${PWD}/petsc-${PETSC_VERSION}
           export PETSC_ARCH=gcc${GCC_V}
           printenv >> $GITHUB_ENV
 
@@ -66,13 +67,13 @@ jobs:
         id: petsc-install
         uses: actions/cache@v2
         with:
-          path: petsc-3.16.0
+          path: petsc-${{ env.PETSC_VERSION }}
-          key: petsc-3.16.0-${{ matrix.os }}-gcc${{ matrix.gcc_v }}-${{ hashFiles('**/petscversion.h') }}
+          key: petsc-${{ env.PETSC_VERSION }}-${{ matrix.os }}-gcc${{ matrix.gcc_v }}-${{ hashFiles('**/petscversion.h') }}
 
       - name: PETSc - Install (Linux)
         if: contains( matrix.os, 'ubuntu')
         run: |
-          cd petsc-3.16.0
+          cd petsc-${PETSC_VERSION}
           ./configure --with-fc=gfortran --with-cc=gcc --with-cxx=g++ \
           --download-mpich --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib \
           --with-mpi-f90module-visibility=0
@@ -81,7 +82,7 @@ jobs:
       - name: PETSc - Install (macOS)
         if: contains( matrix.os, 'macos')
         run: |
-          cd petsc-3.16.0
+          cd petsc-${PETSC_VERSION}
           ./configure --with-fc=gfortran-${GCC_V} --with-cc=gcc-${GCC_V} --with-cxx=g++-${GCC_V} \
           --download-openmpi --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib
           make all
@@ -132,17 +133,17 @@ jobs:
         uses: actions/cache@v2
         with:
           path: download
-          key: petsc-3.16.0.tar.gz
+          key: petsc-${{ env.PETSC_VERSION }}.tar.gz
 
       - name: PETSc - Download
         if: steps.petsc-download.outputs.cache-hit != 'true'
         run: |
-          wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.16.0.tar.gz -P download
+          wget -q https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-${PETSC_VERSION}.tar.gz -P download
 
       - name: PETSc - Prepare
         run: |
-          tar -xf download/petsc-3.16.0.tar.gz -C .
+          tar -xf download/petsc-${PETSC_VERSION}.tar.gz -C .
-          export PETSC_DIR=${PWD}/petsc-3.16.0
+          export PETSC_DIR=${PWD}/petsc-${PETSC_VERSION}
           export PETSC_ARCH=intel-${INTEL_V}
           printenv >> $GITHUB_ENV
 
@@ -150,13 +151,13 @@ jobs:
         id: petsc-install
         uses: actions/cache@v2
         with:
-          path: petsc-3.16.0
+          path: petsc-${{ env.PETSC_VERSION }}
-          key: petsc-3.16.0-intel-${{ matrix.intel_v }}-${{ hashFiles('**/petscversion.h') }}
+          key: petsc-${{ env.PETSC_VERSION }}-intel-${{ matrix.intel_v }}-${{ hashFiles('**/petscversion.h') }}
 
       - name: PETSc - Install (classic)
         if: contains( matrix.intel_v, 'classic')
         run: |
-          cd petsc-3.16.0
+          cd petsc-${PETSC_VERSION}
           ./configure --with-fc=mpiifort --with-cc=mpiicc --with-cxx=mpiicpc \
           --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib
           make all
@@ -164,7 +165,7 @@ jobs:
       - name: PETSc - Install (LLVM)
         if: contains( matrix.intel_v, 'llvm')
         run: |
-          cd petsc-3.16.0
+          cd petsc-${PETSC_VERSION}
           ./configure --with-fc=mpiifort --with-cc="mpiicc -cc=icx" --with-cxx="mpiicpc -cxx=icpx" \
           --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib
           make all

.github/workflows/Python.yml

@@ -43,12 +43,13 @@ jobs:
           pip install pytest
 
       - name: Install dependencies
+        # https://github.com/actions/virtual-environments/issues/4790
         run: >
           sudo apt-get update &&
-          sudo apt-get install python3-pip python3-pytest python3-pandas python3-scipy
+          sudo apt-get remove mysql* &&
-          python3-h5py python3-vtk7 python3-matplotlib python3-yaml -y
+          sudo apt-get install python3-pandas python3-scipy python3-h5py python3-vtk7 python3-matplotlib python3-yaml -y
 
       - name: Run unit tests
         run: |
           export PYTHONPATH=${PWD}/python
-          COLUMNS=256 python -m pytest python
+          COLUMNS=256 pytest python

.gitlab-ci.yml

@@ -36,14 +36,17 @@ variables:
   # Names of module files to load
   # ===============================================================================================
   # ++++++++++++ Compiler +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  COMPILER_INTEL: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
+  COMPILER_GNU:       "Compiler/GNU/10"
-  COMPILER_GNU:   "Compiler/GNU/10"
+  COMPILER_INTELLLVM: "Compiler/oneAPI/2022.0.1 Libraries/IMKL/2022.0.1"
+  COMPILER_INTEL: "Compiler/Intel/2022.0.1 Libraries/IMKL/2022.0.1"
   # ++++++++++++ MPI ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  MPI_INTEL:      "MPI/Intel/19.1.2/IntelMPI/2019"
   MPI_GNU:        "MPI/GNU/10/OpenMPI/4.1.1"
+  MPI_INTELLLVM:  "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0"
+  MPI_INTEL:      "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
   # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  PETSC_INTEL:    "Libraries/PETSc/3.16.1/Intel-19.1.2-IntelMPI-2019"
+  PETSC_GNU:       "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
-  PETSC_GNU:      "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
+  PETSC_INTELLLVM: "Libraries/PETSc/3.16.2/oneAPI-2022.0.1-IntelMPI-2021.5.0"
+  PETSC_INTEL:    "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0"
   # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   MSC:            "FEM/MSC/2021.3.1"
   IntelMarc:      "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
@@ -76,20 +79,6 @@ mypy:
 
 
 ###################################################################################################
-test_grid_Intel:
-  stage: compile
-  script:
-    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
-    - cd PRIVATE/testing/pytest
-    - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel
-
-test_mesh_Intel:
-  stage: compile
-  script:
-    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
-    - cd PRIVATE/testing/pytest
-    - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel
-
 test_grid_GNU:
   stage: compile
   script:
@@ -104,6 +93,27 @@ test_mesh_GNU:
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU
 
+test_mesh_IntelLLVM:
+  stage: compile
+  script:
+    - module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM}
+    - cd PRIVATE/testing/pytest
+    - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_IntelLLVM
+
+test_grid_Intel:
+  stage: compile
+  script:
+    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
+    - cd PRIVATE/testing/pytest
+    - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel
+
+test_mesh_Intel:
+  stage: compile
+  script:
+    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
+    - cd PRIVATE/testing/pytest
+    - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel
+
 test_Marc:
   stage: compile
   script:
@@ -208,7 +218,7 @@ update_revision:
     - git checkout master
     - git pull
     - git merge $UPDATEDREV -s recursive -X ours                                                    # conflicts occur only for inconsistent state
-    - git push origin master                                                                        # master is now tested version and has updated VERSION file
+    - git push
     - git checkout development
     - git pull
     - git merge master -s recursive -X ours -m "[skip ci] Merge branch 'master' into development"   # only possible conflict is in VERSION file

CMakeLists.txt

@@ -82,6 +82,8 @@ if   (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
     include(Compiler-Intel)
 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
     include(Compiler-GNU)
+elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "IntelLLVM")
+    include(Compiler-IntelLLVM)
 else()
     message(FATAL_ERROR "Compiler type(CMAKE_Fortran_COMPILER_ID) not recognized")
 endif()

PRIVATE

@@ -1 +1 @@
-Subproject commit fb4dd2478c00eedf2b5e977b4c2c00dfab523632
+Subproject commit e2cf0c5513b575f2775d2c1b9bfba8c22cd40715

cmake/Compiler-Intel.cmake

@@ -106,8 +106,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0")
 #set (DEBUG_FLAGS "${DEBUG_FLAGS},stderrors")
 #   ... warnings about Fortran standard violations are changed to errors
 
-set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all")
+#set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all")
 # generate debug information for parameters
+# Disabled due to ICE when compiling phase_damage.f90 (not understandable, there is no parameter in there)
 
 # Additional options
 # -heap-arrays:            Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits

cmake/Compiler-IntelLLVM.cmake

@@ -0,0 +1,121 @@
+###################################################################################################
+# Intel Compiler
+###################################################################################################
+if (CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 18.0)
+  message (FATAL_ERROR "Intel Compiler version: ${CMAKE_Fortran_COMPILER_VERSION} not supported")
+endif ()
+
+if (OPENMP)
+  set (OPENMP_FLAGS "-qopenmp")
+endif ()
+
+if (OPTIMIZATION STREQUAL "OFF")
+  set (OPTIMIZATION_FLAGS    "-O0")
+elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
+  set (OPTIMIZATION_FLAGS    "-O2")
+elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
+  set (OPTIMIZATION_FLAGS    "-ipo -O3 -fp-model fast=2 -xHost")
+  # -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost"
+endif ()
+
+# -assume std_mod_proc_name (included in -standard-semantics) causes problems if other modules
+# (PETSc, HDF5) are not compiled with this option (https://software.intel.com/en-us/forums/intel-fortran-compiler-for-linux-and-mac-os-x/topic/62172)
+set (STANDARD_CHECK "-stand f18 -assume nostd_mod_proc_name")
+set (LINKER_FLAGS   "${LINKER_FLAGS} -shared-intel")
+# Link against shared Intel libraries instead of static ones
+
+#------------------------------------------------------------------------------------------------
+# Fine tuning compilation options
+set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp")
+# preprocessor
+
+set (COMPILE_FLAGS "${COMPILE_FLAGS} -ftz")
+# flush underflow to zero, automatically set if -O[1,2,3]
+
+set (COMPILE_FLAGS "${COMPILE_FLAGS} -diag-disable")
+# disables warnings ...
+set (COMPILE_FLAGS "${COMPILE_FLAGS} 5268")
+#   ... the text exceeds right hand column allowed on the line (we have only comments there)
+set (COMPILE_FLAGS "${COMPILE_FLAGS},7624")
+#   ... about deprecated forall (has nice syntax and most likely a performance advantage)
+
+set (COMPILE_FLAGS "${COMPILE_FLAGS} -warn")
+# enables warnings ...
+set (COMPILE_FLAGS "${COMPILE_FLAGS} declarations")
+#   ... any undeclared names (alternative name: -implicitnone)
+set (COMPILE_FLAGS "${COMPILE_FLAGS},general")
+#   ... warning messages and informational messages are issued by the compiler
+set (COMPILE_FLAGS "${COMPILE_FLAGS},usage")
+#   ... questionable programming practices
+set (COMPILE_FLAGS "${COMPILE_FLAGS},interfaces")
+#   ... checks the interfaces of all SUBROUTINEs called and FUNCTIONs invoked in your compilation against an external set of interface blocks
+set (COMPILE_FLAGS "${COMPILE_FLAGS},ignore_loc")
+#   ... %LOC is stripped from an actual argument
+set (COMPILE_FLAGS "${COMPILE_FLAGS},alignments")
+#   ... data that is not naturally aligned
+set (COMPILE_FLAGS "${COMPILE_FLAGS},unused")
+#   ... declared variables that are never used
+
+# Additional options
+#  -warn:                    enables warnings, where
+#     truncated_source:        Determines whether warnings occur when source exceeds the maximum column width in fixed-format files.
+#                               (too many warnings because we have comments beyond character 132)
+#     uncalled:                Determines whether warnings occur when a statement function is never called
+#     all:
+#  -name as_is: case sensitive Fortran!
+
+#------------------------------------------------------------------------------------------------
+# Runtime debugging
+set (DEBUG_FLAGS "${DEBUG_FLAGS} -g")
+# Generate symbolic debugging information in the object file
+
+set (DEBUG_FLAGS "${DEBUG_FLAGS} -traceback")
+# Generate extra information in the object file to provide source file traceback information when a severe error occurs at run time
+
+set (DEBUG_FLAGS "${DEBUG_FLAGS} -gen-interfaces")
+# Generate an interface block for each routine. http://software.intel.com/en-us/blogs/2012/01/05/doctor-fortran-gets-explicit-again/
+
+set (DEBUG_FLAGS "${DEBUG_FLAGS} -fp-stack-check")
+# Generate extra code after every function call to ensure that the floating-point (FP) stack is in the expected state
+
+set (DEBUG_FLAGS "${DEBUG_FLAGS} -fp-model strict")
+# Trap uninitalized variables
+
+set (DEBUG_FLAGS "${DEBUG_FLAGS} -check" )
+# Checks at runtime ...
+set (DEBUG_FLAGS "${DEBUG_FLAGS} bounds")
+#   ... if an array index is too small (<1) or too large!
+set (DEBUG_FLAGS "${DEBUG_FLAGS},format")
+#   ... for the data type of an item being formatted for output.
+set (DEBUG_FLAGS "${DEBUG_FLAGS},output_conversion")
+#   ... for the fit of data items within a designated format descriptor field.
+set (DEBUG_FLAGS "${DEBUG_FLAGS},pointers")
+#   ... for certain disassociated or uninitialized pointers or unallocated allocatable objects.
+set (DEBUG_FLAGS "${DEBUG_FLAGS},uninit")
+#   ... for uninitialized variables.
+set (DEBUG_FLAGS "${DEBUG_FLAGS} -ftrapuv")
+#   ... initializes stack local variables to an unusual value to aid error detection
+set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0")
+#   ... capture all floating-point exceptions, sets -ftz automatically
+
+# disable due to compiler bug https://community.intel.com/t5/Intel-Fortran-Compiler/false-positive-stand-f18-and-IEEE-SELECTED-REAL-KIND/m-p/1227336
+#set (DEBUG_FLAGS "${DEBUG_FLAGS} -warn")
+# enables warnings ...
+#set (DEBUG_FLAGS "${DEBUG_FLAGS} errors")
+#   ... warnings are changed to errors
+#set (DEBUG_FLAGS "${DEBUG_FLAGS},stderrors")
+#   ... warnings about Fortran standard violations are changed to errors
+
+set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all")
+# generate debug information for parameters
+
+# Additional options
+# -heap-arrays:            Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
+# -check:                  Checks at runtime, where
+#    arg_temp_created:       will cause a lot of warnings because we create a bunch of temporary arrays (performance?)
+#    stack:
+
+#------------------------------------------------------------------------------------------------
+#  precision settings
+set (PRECISION_FLAGS "${PRECISION_FLAGS} -real-size 64")
+# set precision for standard real to 32 | 64 | 128 (= 4 | 8 | 16 bytes, type pReal is always 8 bytes)

examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml

@@ -1,17 +1,12 @@
 # Tasan et.al. 2015 Acta Materalia
 # Tasan et.al. 2015 International Journal of Plasticity
 # Diehl et.al. 2015 Meccanica
-Martensite:
+N_sl: [12, 12]
-  lattice: cI
+a_sl: 2.0
-  mechanical:
+dot_gamma_0_sl: 0.001
-    elastic: {C_11: 417.4e+9, C_12: 242.4e+9, C_44: 211.1e+9, type: Hooke}
+h_0_sl-sl: 563.0e+9
-    plastic:
+h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4]
-      N_sl: [12, 12]
+n_sl: 20
-      a_sl: 2.0
+type: phenopowerlaw
-      dot_gamma_0_sl: 0.001
+xi_0_sl: [405.8e+6, 456.7e+6]
-      h_0_sl-sl: 563.0e+9
+xi_inf_sl: [872.9e+6, 971.2e+6]
-      h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4]
-      n_sl: 20
-      type: phenopowerlaw
-      xi_0_sl: [405.8e+6, 456.7e+6]
-      xi_inf_sl: [872.9e+6, 971.2e+6]

examples/config/phase/AISI304.yaml

@@ -0,0 +1,6 @@
+references:
+  - H.M. Ledbetter
+    physica status solidi (a) 85(1):89-96, 1984
+    https://doi.org/10.1002/pssa.2210850111
+lattice: cF
+rho: 7937.0

examples/config/phase/Ag.yaml

@@ -0,0 +1,4 @@
+references:
+  - https://en.wikipedia.org/wiki/Silver
+lattice: cF
+rho: 10490.0

examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml

@@ -0,0 +1,7 @@
+type: thermalexpansion
+references:
+  - R.H. Bogaard et al.
+    Thermochimica Acta 218:373-393, 1993
+    https://doi.org/10.1016/0040-6031(93)80437-F
+A_11: 15.0e-6
+T_ref: 300.0

examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml

@@ -0,0 +1,8 @@
+type: Hooke
+references:
+  - H.M. Ledbetter
+    physica status solidi (a) 85(1):89-96, 1984
+    https://doi.org/10.1002/pssa.2210850111
+C_11: 204.6e+9
+C_12: 137.7e+9
+C_44: 126.2e+9

examples/config/phase/mechanical/elastic/Hooke_Ag.yaml

@@ -0,0 +1,22 @@
+type: Hooke
+references:
+  - J.R. Neighbours and G.A. Alers,
+    Physical Review 111:707-712, 1958,
+    https://doi.org/10.1103/PhysRev.111.707
+  - Y.A. Chang and L. Himmel,
+    Journal of Applied Physics 37:3567-3572, 1966,
+    https://doi.org/10.1063/1.1708903
+
+T_ref: 300
+
+C_11: 122.9e+9
+C_11,T: -313.5e+5
+C_11,T^2: -107.3e+2
+
+C_12: 91.55e+9
+C_12,T: -164.1e+5
+C_12,T^2: -681.6e+1
+
+C_44: 42.63e+9
+C_44,T: -180.5e+5
+C_44,T^2: -353.8e+1

examples/config/phase/mechanical/elastic/Hooke_Al.yaml

@@ -1,8 +1,22 @@
 type: Hooke
 references:
-  - J. Vallin et al.,
+  - G.N. Kamm and G.A. Alers,
-    Journal of Applied Physics 35(6):1825-1826, 1964,
+    Journal of Applied Physics 35:327-330, 1964,
-    https://doi.org/10.1063/1.1713749
+    https://doi.org/10.1063/1.1713309
-C_11: 107.3e+9
+  - D. Gerlich and E.S. Fisher,
-C_12: 60.8e+9
+    Journal of Physics and Chemistry of Solids 30:1197-1205, 1969
-C_44: 28.3e+9
+    https://doi.org/10.1016/0022-3697(69)90377-1
+
+T_ref: 300
+
+C_11: 106.1e+9
+C_11,T: -359.3e+5
+C_11,T^2: -152.7e+2
+
+C_12: 57.83e+9
+C_12,T: -781.6e+4
+C_12,T^2: -551.3e+1
+
+C_44: 24.31e+9
+C_44,T: -142.9e+5
+C_44,T^2: -404.6e+1

examples/config/phase/mechanical/elastic/Hooke_Fe.yaml

@@ -1,9 +1,19 @@
 type: Hooke
 references:
-  - J.P. Hirth and J. Lothe,
+  - D.J. Dever,
-    Theory of Dislocations, 1982,
+    Journal of Applied Physics 43(8):3293-3301, 1972,
-    John Wiley & Sons,
+    https://doi.org/10.1063/1.1661710
-    page 837
+
-C_11: 242.e9
+T_ref: 300
-C_12: 146.5e+9
+
-C_44: 112.e9
+C_11: 231.7e+9
+C_11,T: -47.6e+6
+C_11,T^2: -54.4e+3
+
+C_12: 135.8e+9
+C_12,T: -12.9e+6
+C_12,T^2: -7.3e+3
+
+C_44: 116.8e+9
+C_44,T: -19.4e+6
+C_44,T^2: -2.5e+3

examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml

@@ -0,0 +1,8 @@
+type: Hooke
+references:
+  - S.A. Kim and W.L. Johnson,
+    Materials Science & Engineering A 452-453:633-639, 2007,
+    https://doi.org/10.1016/j.msea.2006.11.147
+C_11: 268.1e+9
+C_12: 111.2e+9
+C_44: 79.06e+9

examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml

@@ -4,7 +4,8 @@ references:
     International Journal of Plasticity 134:102779, 2020,
     https://doi.org/10.1016/j.ijplas.2020.102779
   - K. Sedighiani et al.,
-    Mechanics of Materials, submitted
+    Mechanics of Materials, 164:104117, 2022,
+    https://doi.org/10.1016/j.mechmat.2021.104117
 output: [rho_dip, rho_mob]
 N_sl: [12, 12]
 b_sl: [2.49e-10, 2.49e-10]

examples/config/phase/thermal/AISI304.yaml

@@ -0,0 +1,9 @@
+references:
+  - B.F. Blackwell et al.
+    Proceedings of 34th National Heat Transfer Conference 2000
+    https://www.osti.gov/servlets/purl/760791
+  - R.H. Bogaard et al.
+    Thermochimica Acta 218:373-393, 1993
+    https://doi.org/10.1016/0040-6031(93)80437-F
+C_p: 470.0
+K_11: 14.34

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha5-244-gb57351045
+v3.0.0-alpha5-358-g81a7c32a5

python/damask/_colormap.py

@@ -6,6 +6,11 @@ from pathlib import Path
 from typing import Sequence, Union, TextIO
 
 import numpy as np
+try:
+    from numpy.typing import ArrayLike
+except ImportError:
+    ArrayLike = Union[np.ndarray,Sequence[float]] # type: ignore
+import scipy.interpolate as interp
 import matplotlib as mpl
 if os.name == 'posix' and 'DISPLAY' not in os.environ:
     mpl.use('Agg')
@@ -41,7 +46,7 @@ class Colormap(mpl.colors.ListedColormap):
     https://doi.org/10.1016/j.ijplas.2012.09.012
 
     Matplotlib colormaps overview
-    https://matplotlib.org/tutorials/colors/colormaps.html
+    https://matplotlib.org/stable/tutorials/colors/colormaps.html
 
     """
 
@@ -77,8 +82,8 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def from_range(low: Sequence[float],
+    def from_range(low: ArrayLike,
-                   high: Sequence[float],
+                   high: ArrayLike,
                    name: str = 'DAMASK colormap',
                    N: int = 256,
                    model: str = 'rgb') -> 'Colormap':
@@ -128,7 +133,7 @@ class Colormap(mpl.colors.ListedColormap):
         if model.lower() not in toMsh:
             raise ValueError(f'Invalid color model: {model}.')
 
-        low_high = np.vstack((low,high))
+        low_high = np.vstack((low,high)).astype(float)
         out_of_bounds = np.bool_(False)
 
         if   model.lower() == 'rgb':
@@ -141,7 +146,7 @@ class Colormap(mpl.colors.ListedColormap):
             out_of_bounds = np.any(low_high[:,0]<0)
 
         if out_of_bounds:
-            raise ValueError(f'{model.upper()} colors {low} | {high} are out of bounds.')
+            raise ValueError(f'{model.upper()} colors {low_high[0]} | {low_high[1]} are out of bounds.')
 
         low_,high_ = map(toMsh[model.lower()],low_high)
         msh = map(functools.partial(Colormap._interpolate_msh,low=low_,high=high_),np.linspace(0,1,N))
@@ -191,19 +196,50 @@ class Colormap(mpl.colors.ListedColormap):
             return Colormap.from_range(definition['low'],definition['high'],name,N)
 
 
+    def at(self,
+           fraction : Union[float,Sequence[float]]) -> np.ndarray:
+        """
+        Interpolate color at fraction.
+
+        Parameters
+        ----------
+        fraction : float or sequence of float
+            Fractional coordinate(s) to evaluate Colormap at.
+
+        Returns
+        -------
+        color : np.ndarray, shape(...,4)
+            RGBA values of interpolated color(s).
+
+        Examples
+        --------
+        >>> import damask
+        >>> cmap = damask.Colormap.from_predefined('gray')
+        >>> cmap.at(0.5)
+        array([0.5, 0.5, 0.5, 1. ])
+        >>> 'rgb({},{},{})'.format(*cmap.at(0.5))
+        'rgb(0.5,0.5,0.5)'
+
+        """
+        return interp.interp1d(np.linspace(0,1,self.N),
+                               self.colors,
+                               axis=0,
+                               assume_sorted=True)(fraction)
+
+
     def shade(self,
               field: np.ndarray,
-              bounds: Sequence[float] = None,
+              bounds: ArrayLike = None,
               gap: float = None) -> Image:
         """
         Generate PIL image of 2D field using colormap.
 
         Parameters
         ----------
-        field : numpy.array, shape (:,:)
+        field : numpy.ndarray, shape (:,:)
             Data to be shaded.
         bounds : sequence of float, len (2), optional
-            Value range (low,high) spanned by colormap.
+            Value range (left,right) spanned by colormap.
         gap : field.dtype, optional
             Transparent value. NaN will always be rendered transparent.
 
@@ -213,21 +249,20 @@ class Colormap(mpl.colors.ListedColormap):
             RGBA image of shaded data.
 
         """
-        N = len(self.colors)
         mask = np.logical_not(np.isnan(field) if gap is None else \
                np.logical_or (np.isnan(field), field == gap))                                       # mask NaN (and gap if present)
 
-        lo,hi = (field[mask].min(),field[mask].max()) if bounds is None else \
+        l,r = (field[mask].min(),field[mask].max()) if bounds is None else \
-                (min(bounds[:2]),max(bounds[:2]))
+              np.array(bounds,float)[:2]
 
-        delta,avg = hi-lo,0.5*(hi+lo)
+        delta,avg = r-l,0.5*abs(r+l)
 
-        if delta * 1e8 <= avg:                                                                      # delta is similar to numerical noise
+        if abs(delta) * 1e8 <= avg:                                                                 # delta is similar to numerical noise
-            hi,lo = hi+0.5*avg,lo-0.5*avg                                                           # extend range to have actual data centered within
+            l,r = l-0.5*avg*np.sign(delta),r+0.5*avg*np.sign(delta),                                # extend range to have actual data centered within
 
         return Image.fromarray(
             (np.dstack((
-                        self.colors[(np.round(np.clip((field-lo)/(hi-lo),0.0,1.0)*(N-1))).astype(np.uint16),:3],
+                        self.colors[(np.round(np.clip((field-l)/delta,0.0,1.0)*(self.N-1))).astype(np.uint16),:3],
                         mask.astype(float)
                        )
                       )*255
@@ -343,7 +378,7 @@ class Colormap(mpl.colors.ListedColormap):
         # ToDo: test in GOM
         GOM_str = '1 1 {name} 9 {name} '.format(name=self.name.replace(" ","_")) \
                 +  '0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 ' \
-                + f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {len(self.colors)}' \
+                + f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {self.N}' \
                 + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(int))]) \
                 + '\n'
 
@@ -581,7 +616,7 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
+    def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
         """
         CIE Lab to CIE Xyz.
 
@@ -589,6 +624,8 @@ class Colormap(mpl.colors.ListedColormap):
         ----------
         lab : numpy.ndarray, shape (3)
             CIE lab values.
+        ref_white : numpy.ndarray, shape (3)
+            Reference white, default value is the standard 2° observer for D65.
 
         Returns
         -------
@@ -607,10 +644,10 @@ class Colormap(mpl.colors.ListedColormap):
                          f_x**3.                if f_x**3. > _EPS     else (116.*f_x-16.)/_KAPPA,
                          ((lab[0]+16.)/116.)**3 if lab[0]>_KAPPA*_EPS else lab[0]/_KAPPA,
                          f_z**3.                if f_z**3. > _EPS     else (116.*f_z-16.)/_KAPPA
-                        ])*(ref_white if ref_white is not None else _REF_WHITE)
+                        ])*ref_white
 
     @staticmethod
-    def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
+    def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
         """
         CIE Xyz to CIE Lab.
 
@@ -618,6 +655,8 @@ class Colormap(mpl.colors.ListedColormap):
         ----------
         xyz : numpy.ndarray, shape (3)
             CIE Xyz values.
+        ref_white : numpy.ndarray, shape (3)
+            Reference white, default value is the standard 2° observer for D65.
 
         Returns
         -------
@@ -629,7 +668,6 @@ class Colormap(mpl.colors.ListedColormap):
         http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
 
         """
-        ref_white = ref_white if ref_white is not None else _REF_WHITE
         f = np.where(xyz/ref_white > _EPS,(xyz/ref_white)**(1./3.),(_KAPPA*xyz/ref_white+16.)/116.)
 
         return np.array([

python/damask/_crystal.py

@@ -1,3 +1,5 @@
+from typing import Union, Dict, List, Tuple
+
 import numpy as np
 
 from . import util
@@ -177,7 +179,7 @@ class Crystal():
 
 
     @property
-    def standard_triangle(self):
+    def standard_triangle(self) -> Union[Dict[str, np.ndarray], None]:
         """
         Corners of the standard triangle.
 
@@ -238,7 +240,7 @@ class Crystal():
                                                [-1., 0., 0.],
                                                [ 0., 1., 0.] ]),
                         }}
-        return _basis.get(self.family,None)
+        return _basis.get(self.family, None)
 
 
     @property
@@ -256,7 +258,7 @@ class Crystal():
 
 
     @property
-    def basis_real(self):
+    def basis_real(self) -> np.ndarray:
         """
         Return orthogonal real space crystal basis.
 
@@ -280,7 +282,7 @@ class Crystal():
 
 
     @property
-    def basis_reciprocal(self):
+    def basis_reciprocal(self) -> np.ndarray:
         """Return reciprocal (dual) crystal basis."""
         return np.linalg.inv(self.basis_real.T)
 
@@ -312,7 +314,7 @@ class Crystal():
                         + _lattice_points.get(self.lattice if self.lattice == 'hP' else \
                                               self.lattice[-1],None),dtype=float)
 
-    def to_lattice(self,*,direction=None,plane=None):
+    def to_lattice(self, *, direction: np.ndarray = None, plane: np.ndarray = None) -> np.ndarray:
         """
         Calculate lattice vector corresponding to crystal frame direction or plane normal.
 
@@ -336,7 +338,7 @@ class Crystal():
         return np.einsum('il,...l',basis,axis)
 
 
-    def to_frame(self,*,uvw=None,hkl=None):
+    def to_frame(self, *, uvw: np.ndarray = None, hkl: np.ndarray = None) -> np.ndarray:
         """
         Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl).
 
@@ -359,7 +361,7 @@ class Crystal():
         return np.einsum('il,...l',basis,axis)
 
 
-    def kinematics(self,mode):
+    def kinematics(self, mode: str) -> Dict[str, List[np.ndarray]]:
         """
         Return crystal kinematics systems.
 
@@ -551,7 +553,7 @@ class Crystal():
                     'plane':    [m[:,3:6] for m in master]}
 
 
-    def relation_operations(self,model):
+    def relation_operations(self, model: str) -> Tuple[str, Rotation]:
         """
         Crystallographic orientation relationships for phase transformations.
 

python/damask/_grid.py

@@ -197,7 +197,7 @@ class Grid:
             Grid-based geometry from file.
 
         """
-        warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.1.0', DeprecationWarning,2)
+        warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.0.0', DeprecationWarning,2)
         try:
             f = open(fname)
         except TypeError:
@@ -629,7 +629,7 @@ class Grid:
             Compress geometry with 'x of y' and 'a to b'.
 
         """
-        warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.1.0', DeprecationWarning,2)
+        warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.0.0', DeprecationWarning,2)
         header =  [f'{len(self.comments)+4} header'] + self.comments \
                 + ['grid   a {} b {} c {}'.format(*self.cells),
                    'size   x {} y {} z {}'.format(*self.size),

python/damask/_orientation.py

@@ -13,7 +13,7 @@ _parameter_doc = \
        """
         family : {'triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic'}, optional.
             Name of the crystal family.
-            Will be infered if 'lattice' is given.
+            Family will be inferred if 'lattice' is given.
         lattice : {'aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF'}, optional.
             Name of the Bravais lattice in Pearson notation.
         a : float, optional

python/damask/_result.py

@@ -4,6 +4,7 @@ import fnmatch
 import os
 import copy
 import datetime
+import warnings
 import xml.etree.ElementTree as ET
 import xml.dom.minidom
 from pathlib import Path
@@ -27,6 +28,20 @@ h5py3 = h5py.__version__[0] == '3'
 
 chunk_size = 1024**2//8                                                                             # for compression in HDF5
 
+def _view_transition(what,datasets,increments,times,phases,homogenizations,fields):
+    if (datasets is not None and what is None) or (what is not None and datasets is None):
+        raise ValueError('"what" and "datasets" need to be used as a pair')
+    if datasets is not None or what is not None:
+        warnings.warn('Arguments "what" and "datasets" will be removed in DAMASK v3.0.0-alpha7', DeprecationWarning,2)
+        return what,datasets
+    if sum(1 for _ in filter(None.__ne__, [increments,times,phases,homogenizations,fields])) > 1:
+        raise ValueError('Only one out of "increments", "times", "phases", "homogenizations", and "fields" can be used')
+    else:
+        if increments is not None: return "increments", increments
+        if times is not None: return "times", times
+        if phases is not None: return "phases", phases
+        if homogenizations is not None: return "homogenizations", homogenizations
+        if fields is not None: return "fields", fields
 
 def _read(dataset):
     """Read a dataset and its metadata into a numpy.ndarray."""
@@ -79,7 +94,7 @@ class Result:
     >>> r.add_Cauchy()
     >>> r.add_equivalent_Mises('sigma')
     >>> r.export_VTK()
-    >>> r_last = r.view('increments',-1)
+    >>> r_last = r.view(increments=-1)
     >>> sigma_vM_last = r_last.get('sigma_vM')
 
     """
@@ -141,7 +156,7 @@ class Result:
 
         self.fname = Path(fname).absolute()
 
-        self._allow_modification = False
+        self._protected = True
 
 
     def __copy__(self):
@@ -155,10 +170,10 @@ class Result:
         """Show summary of file content."""
         visible_increments = self.visible['increments']
 
-        first = self.view('increments',visible_increments[0:1]).list_data()
+        first = self.view(increments=visible_increments[0:1]).list_data()
 
         last  = '' if len(visible_increments) < 2 else \
-                self.view('increments',visible_increments[-1:]).list_data()
+                self.view(increments=visible_increments[-1:]).list_data()
 
         in_between = '' if len(visible_increments) < 3 else \
                      ''.join([f'\n{inc}\n  ...\n' for inc in visible_increments[1:-1]])
@@ -231,36 +246,6 @@ class Result:
         return dup
 
 
-    def modification_enable(self):
-        """
-        Allow modification of existing data.
-
-        Returns
-        -------
-        modified_view : damask.Result
-            View without write-protection of existing data.
-
-        """
-        print(util.warn('Warning: Modification of existing datasets allowed!'))
-        dup = self.copy()
-        dup._allow_modification = True
-        return dup
-
-    def modification_disable(self):
-        """
-        Prevent modification of existing data (default case).
-
-        Returns
-        -------
-        modified_view : damask.Result
-            View with write-protection of existing data.
-
-        """
-        dup = self.copy()
-        dup._allow_modification = False
-        return dup
-
-
     def increments_in_range(self,start,end):
         """
         Get all increments within a given range.
@@ -285,7 +270,6 @@ class Result:
                 selected.append(self.increments[i])
         return selected
 
-
     def times_in_range(self,start,end):
         """
         Get all increments within a given time range.
@@ -310,17 +294,38 @@ class Result:
         return selected
 
 
-    def view(self,what,datasets):
+    def view(self,what=None,datasets=None,*,
+                  increments=None,
+                  times=None,
+                  phases=None,
+                  homogenizations=None,
+                  fields=None,
+                  protected=None):
         """
         Set view.
 
+        Wildcard matching with '?' and '*' is supported.
+        True is equivalent to '*', False is equivalent to [].
+
         Parameters
         ----------
         what : {'increments', 'times', 'phases', 'homogenizations', 'fields'}
-            Attribute to change.
+            Attribute to change. DEPRECATED.
         datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool
-            Name of datasets; supports '?' and '*' wildcards.
+            Name of datasets; supports '?' and '*' wildcards. DEPRECATED.
             True is equivalent to '*', False is equivalent to [].
+        increments: (list of) int, (list of) str, or bool, optional.
+            Number(s) of increments to select.
+        times: (list of) float, (list of) str, or bool, optional.
+            Simulation time(s) of increments to select.
+        phases: (list of) str, or bool, optional.
+            Name(s) of phases to select.
+        homogenizations: (list of) str, or bool, optional.
+            Name(s) of homogenizations to select.
+        fields: (list of) str, or bool, optional.
+            Name(s) of fields to select.
+        protected: bool, optional.
+            Protection status of existing data.
 
         Returns
         -------
@@ -333,29 +338,61 @@ class Result:
 
         >>> import damask
         >>> r = damask.Result('my_file.hdf5')
-        >>> r_first = r.view('increment',0)
+        >>> r_first = r.view(increment=0)
 
         Get a view that shows all results between simulation times of 10 to 40:
 
         >>> import damask
         >>> r = damask.Result('my_file.hdf5')
-        >>> r_t10to40 = r.view('times',r.times_in_range(10.0,40.0))
+        >>> r_t10to40 = r.view(times=r.times_in_range(10.0,40.0))
 
         """
-        return self._manage_view('set',what,datasets)
+        v = _view_transition(what,datasets,increments,times,phases,homogenizations,fields)
+        if protected is not None:
+            if v is None:
+                dup = self.copy()
+            else:
+                what_,datasets_ = v
+                dup = self._manage_view('set',what_,datasets_)
+            if not protected:
+                print(util.warn('Warning: Modification of existing datasets allowed!'))
+            dup._protected = protected
+        else:
+            what_,datasets_ = v
+            dup = self._manage_view('set',what_,datasets_)
+
+        return dup
 
 
-    def view_more(self,what,datasets):
+    def view_more(self,what=None,datasets=None,*,
+                  increments=None,
+                  times=None,
+                  phases=None,
+                  homogenizations=None,
+                  fields=None):
         """
         Add to view.
 
+        Wildcard matching with '?' and '*' is supported.
+        True is equivalent to '*', False is equivalent to [].
+
         Parameters
         ----------
         what : {'increments', 'times', 'phases', 'homogenizations', 'fields'}
-            Attribute to change.
+            Attribute to change. DEPRECATED.
         datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool
-            Name of datasets; supports '?' and '*' wildcards.
+            Name of datasets; supports '?' and '*' wildcards. DEPRECATED.
             True is equivalent to '*', False is equivalent to [].
+        increments: (list of) int, (list of) str, or bool, optional.
+            Number(s) of increments to select.
+        times: (list of) float, (list of) str, or bool, optional.
+            Simulation time(s) of increments to select.
+        phases: (list of) str, or bool, optional.
+            Name(s) of phases to select.
+        homogenizations: (list of) str, or bool, optional.
+            Name(s) of homogenizations to select.
+        fields: (list of) str, or bool, optional.
+            Name(s) of fields to select.
 
         Returns
         -------
@@ -367,25 +404,44 @@ class Result:
         Get a view that shows only results from first and last increment:
 
         >>> import damask
-        >>> r_empty = damask.Result('my_file.hdf5').view('increments',False)
+        >>> r_empty = damask.Result('my_file.hdf5').view(increments=False)
-        >>> r_first = r_empty.view_more('increments',0)
+        >>> r_first = r_empty.view_more(increments=0)
-        >>> r_first_and_last = r.first.view_more('increments',-1)
+        >>> r_first_and_last = r.first.view_more(increments=-1)
 
         """
-        return self._manage_view('add',what,datasets)
+        what_, datasets_ = _view_transition(what,datasets,increments,times,phases,homogenizations,fields)
+        return self._manage_view('add',what_,datasets_)
 
 
-    def view_less(self,what,datasets):
+    def view_less(self,what=None,datasets=None,*,
+                  increments=None,
+                  times=None,
+                  phases=None,
+                  homogenizations=None,
+                  fields=None):
         """
         Remove from view.
 
+        Wildcard matching with '?' and '*' is supported.
+        True is equivalent to '*', False is equivalent to [].
+
         Parameters
         ----------
         what : {'increments', 'times', 'phases', 'homogenizations', 'fields'}
-            Attribute to change.
+            Attribute to change. DEPRECATED.
         datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool
-            Name of datasets; supports '?' and '*' wildcards.
+            Name of datasets; supports '?' and '*' wildcards. DEPRECATED.
             True is equivalent to '*', False is equivalent to [].
+        increments: (list of) int, (list of) str, or bool, optional.
+            Number(s) of increments to select.
+        times: (list of) float, (list of) str, or bool, optional.
+            Simulation time(s) of increments to select.
+        phases: (list of) str, or bool, optional.
+            Name(s) of phases to select.
+        homogenizations: (list of) str, or bool, optional.
+            Name(s) of homogenizations to select.
+        fields: (list of) str, or bool, optional.
+            Name(s) of fields to select.
 
         Returns
         -------
@@ -398,10 +454,11 @@ class Result:
 
         >>> import damask
         >>> r_all = damask.Result('my_file.hdf5')
-        >>> r_deformed = r_all.view_less('increments',0)
+        >>> r_deformed = r_all.view_less(increments=0)
 
         """
-        return self._manage_view('del',what,datasets)
+        what_, datasets_ = _view_transition(what,datasets,increments,times,phases,homogenizations,fields)
+        return self._manage_view('del',what_,datasets_)
 
 
     def rename(self,name_src,name_dst):
@@ -424,11 +481,11 @@ class Result:
 
         >>> import damask
         >>> r = damask.Result('my_file.hdf5')
-        >>> r_unprotected = r.modification_enable()
+        >>> r_unprotected = r.view(protected=False)
         >>> r_unprotected.rename('F','def_grad')
 
         """
-        if not self._allow_modification:
+        if self._protected:
             raise PermissionError('Renaming datasets not permitted')
 
         with h5py.File(self.fname,'a') as f:
@@ -463,11 +520,11 @@ class Result:
 
         >>> import damask
         >>> r = damask.Result('my_file.hdf5')
-        >>> r_unprotected = r.modification_enable()
+        >>> r_unprotected = r.view(protected=False)
         >>> r_unprotected.remove('F')
 
         """
-        if not self._allow_modification:
+        if self._protected:
             raise PermissionError('Removing datasets not permitted')
 
         with h5py.File(self.fname,'a') as f:
@@ -1358,7 +1415,7 @@ class Result:
             lock.acquire()
             with h5py.File(self.fname, 'a') as f:
                 try:
-                    if self._allow_modification and '/'.join([group,result['label']]) in f:
+                    if not self._protected and '/'.join([group,result['label']]) in f:
                         dataset = f['/'.join([group,result['label']])]
                         dataset[...] = result['data']
                         dataset.attrs['overwritten'] = True

python/damask/util.py

@@ -18,11 +18,11 @@ from . import version
 __all__=[
          'srepr',
          'emph', 'deemph', 'warn', 'strikeout',
-         'run', 
+         'run',
          'natural_sort',
          'show_progress',
          'scale_to_coprime',
-         'project_stereographic',
+         'project_equal_angle', 'project_equal_area',
          'hybrid_IA',
          'execution_stamp',
          'shapeshifter', 'shapeblender',
@@ -267,13 +267,13 @@ def scale_to_coprime(v):
     return m
 
 
-def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
+def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
     """
-    Apply stereographic projection to vector.
+    Apply equal-angle projection to vector.
 
     Parameters
     ----------
-    vector : numpy.ndarray of shape (...,3)
+    vector : numpy.ndarray, shape (...,3)
         Vector coordinates to be projected.
     direction : str
         Projection direction 'x', 'y', or 'z'.
@@ -281,32 +281,74 @@ def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
     normalize : bool
         Ensure unit length of input vector. Defaults to True.
     keepdims : bool
-        Maintain three-dimensional output coordinates.
+        Maintain three-dimensional output coordinates. Defaults to False.
-        Default two-dimensional output uses right-handed frame spanned by
+        Two-dimensional output uses right-handed frame spanned by
         the next and next-next axis relative to the projection direction,
         e.g. x-y when projecting along z and z-x when projecting along y.
 
     Returns
     -------
-    coordinates : numpy.ndarray of shape (...,2 | 3)
+    coordinates : numpy.ndarray, shape (...,2 | 3)
         Projected coordinates.
 
     Examples
     --------
     >>> import damask
     >>> import numpy as np
-    >>> project_stereographic(np.ones(3))
+    >>> project_equal_angle(np.ones(3))
         [0.3660254, 0.3660254]
-    >>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True)
+    >>> project_equal_angle(np.ones(3),direction='x',normalize=False,keepdims=True)
         [0, 0.5, 0.5]
-    >>> project_stereographic([0,1,1],direction='y',normalize=True,keepdims=False)
+    >>> project_equal_angle([0,1,1],direction='y',normalize=True,keepdims=False)
         [0.41421356, 0]
 
     """
     shift = 'zyx'.index(direction)
-    v_ = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
+    v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
-                 shift,axis=-1)
+                shift,axis=-1)
-    return np.roll(np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])),np.zeros_like(v_[...,2:3])]),
+    return np.roll(np.block([v[...,:2]/(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
+                   -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
+
+def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
+    """
+    Apply equal-area projection to vector.
+
+    Parameters
+    ----------
+    vector : numpy.ndarray, shape (...,3)
+        Vector coordinates to be projected.
+    direction : str
+        Projection direction 'x', 'y', or 'z'.
+        Defaults to 'z'.
+    normalize : bool
+        Ensure unit length of input vector. Defaults to True.
+    keepdims : bool
+        Maintain three-dimensional output coordinates. Defaults to False.
+        Two-dimensional output uses right-handed frame spanned by
+        the next and next-next axis relative to the projection direction,
+        e.g. x-y when projecting along z and z-x when projecting along y.
+
+    Returns
+    -------
+    coordinates : numpy.ndarray, shape (...,2 | 3)
+        Projected coordinates.
+
+    Examples
+    --------
+    >>> import damask
+    >>> import numpy as np
+    >>> project_equal_area(np.ones(3))
+        [0.45970084, 0.45970084]
+    >>> project_equal_area(np.ones(3),direction='x',normalize=False,keepdims=True)
+        [0.0, 0.70710678, 0.70710678]
+    >>> project_equal_area([0,1,1],direction='y',normalize=True,keepdims=False)
+        [0.5411961, 0.0]
+
+    """
+    shift = 'zyx'.index(direction)
+    v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
+                shift,axis=-1)
+    return np.roll(np.block([v[...,:2]/np.sqrt(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
                    -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
 
 

python/tests/test_Colormap.py

@@ -139,6 +139,16 @@ class TestColormap:
         c += c
         assert (np.allclose(c.colors[:len(c.colors)//2],c.colors[len(c.colors)//2:]))
 
+    @pytest.mark.parametrize('N,cmap,at,result',[
+           (8,'gray',0.5,[0.5,0.5,0.5]),
+           (17,'gray',0.5,[0.5,0.5,0.5]),
+           (17,'gray',[0.5,0.75],[[0.5,0.5,0.5],[0.75,0.75,0.75]]),
+           ])
+    def test_at_value(self, N, cmap, at, result):
+        assert np.allclose(Colormap.from_predefined(cmap,N=N).at(at)[...,:3],
+                           result,
+                           rtol=0.005)
+
     @pytest.mark.parametrize('bounds',[None,[2,10]])
     def test_shade(self,ref_path,update,bounds):
         data = np.add(*np.indices((10, 11)))

python/tests/test_Result.py

@@ -25,7 +25,7 @@ def default(tmp_path,ref_path):
     fname = '12grains6x7x8_tensionY.hdf5'
     shutil.copy(ref_path/fname,tmp_path)
     f = Result(tmp_path/fname)
-    return f.view('times',20.0)
+    return f.view(times=20.0)
 
 @pytest.fixture
 def single_phase(tmp_path,ref_path):
@@ -58,14 +58,14 @@ class TestResult:
 
 
     def test_view_all(self,default):
-        a = default.view('increments',True).get('F')
+        a = default.view(increments=True).get('F')
 
-        assert dict_equal(a,default.view('increments','*').get('F'))
+        assert dict_equal(a,default.view(increments='*').get('F'))
-        assert dict_equal(a,default.view('increments',default.increments_in_range(0,np.iinfo(int).max)).get('F'))
+        assert dict_equal(a,default.view(increments=default.increments_in_range(0,np.iinfo(int).max)).get('F'))
 
-        assert dict_equal(a,default.view('times',True).get('F'))
+        assert dict_equal(a,default.view(times=True).get('F'))
-        assert dict_equal(a,default.view('times','*').get('F'))
+        assert dict_equal(a,default.view(times='*').get('F'))
-        assert dict_equal(a,default.view('times',default.times_in_range(0.0,np.inf)).get('F'))
+        assert dict_equal(a,default.view(times=default.times_in_range(0.0,np.inf)).get('F'))
 
     @pytest.mark.parametrize('what',['increments','times','phases','fields'])                       # ToDo: discuss homogenizations
     def test_view_none(self,default,what):
@@ -314,7 +314,7 @@ class TestResult:
 
     @pytest.mark.parametrize('overwrite',['off','on'])
     def test_add_overwrite(self,default,overwrite):
-        last = default.view('increments',-1)
+        last = default.view(increments=-1)
 
         last.add_stress_Cauchy()
 
@@ -322,9 +322,9 @@ class TestResult:
         created_first = datetime.strptime(created_first,'%Y-%m-%d %H:%M:%S%z')
 
         if overwrite == 'on':
-            last = last.modification_enable()
+            last = last.view(protected=False)
         else:
-            last = last.modification_disable()
+            last = last.view(protected=True)
 
         time.sleep(2.)
         try:
@@ -344,10 +344,10 @@ class TestResult:
     def test_rename(self,default,allowed):
         if allowed == 'on':
             F = default.place('F')
-            default = default.modification_enable()
+            default = default.view(protected=False)
             default.rename('F','new_name')
             assert np.all(F == default.place('new_name'))
-            default = default.modification_disable()
+            default = default.view(protected=True)
 
         with pytest.raises(PermissionError):
             default.rename('P','another_new_name')
@@ -355,7 +355,7 @@ class TestResult:
     @pytest.mark.parametrize('allowed',['off','on'])
     def test_remove(self,default,allowed):
         if allowed == 'on':
-            unsafe = default.modification_enable()
+            unsafe = default.view(protected=False)
             unsafe.remove('F')
             assert unsafe.get('F') is None
         else:
@@ -377,7 +377,7 @@ class TestResult:
     @pytest.mark.parametrize('inc',[4,0],ids=range(2))
     @pytest.mark.xfail(int(vtk.vtkVersion.GetVTKVersion().split('.')[0])<9, reason='missing "Direction" attribute')
     def test_vtk(self,request,tmp_path,ref_path,update,patch_execution_stamp,patch_datetime_now,output,fname,inc):
-        result = Result(ref_path/fname).view('increments',inc)
+        result = Result(ref_path/fname).view(increments=inc)
         os.chdir(tmp_path)
         result.export_VTK(output,parallel=False)
         fname = fname.split('.')[0]+f'_inc{(inc if type(inc) == int else inc[0]):0>2}.vti'
@@ -400,7 +400,7 @@ class TestResult:
         result.export_VTK(output,mode)
 
     def test_marc_coordinates(self,ref_path):
-        result = Result(ref_path/'check_compile_job1.hdf5').view('increments',-1)
+        result = Result(ref_path/'check_compile_job1.hdf5').view(increments=-1)
         c_n = result.coordinates0_node + result.get('u_n')
         c_p = result.coordinates0_point + result.get('u_p')
         assert len(c_n) > len(c_p)
@@ -440,7 +440,7 @@ class TestResult:
         dim_xdmf = reader_xdmf.GetOutput().GetDimensions()
         bounds_xdmf = reader_xdmf.GetOutput().GetBounds()
 
-        single_phase.view('increments',0).export_VTK(parallel=False)
+        single_phase.view(increments=0).export_VTK(parallel=False)
         fname = os.path.splitext(os.path.basename(single_phase.fname))[0]+'_inc00.vti'
         reader_vti = vtk.vtkXMLImageDataReader()
         reader_vti.SetFileName(fname)

python/tests/test_VTK.py

@@ -20,7 +20,7 @@ def default():
     """Simple VTK."""
     cells = np.array([5,6,7],int)
     size  = np.array([.6,1.,.5])
-    return VTK.from_rectilinear_grid(cells,size)
+    return VTK.from_image_data(cells,size)
 
 class TestVTK:
 
@@ -116,7 +116,7 @@ class TestVTK:
 
     def test_add_extension(self,tmp_path,default):
         default.save(tmp_path/'default.txt',parallel=False)
-        assert os.path.isfile(tmp_path/'default.txt.vtr')
+        assert os.path.isfile(tmp_path/'default.txt.vti')
 
 
     def test_invalid_get(self,default):
@@ -160,7 +160,7 @@ class TestVTK:
     def test_comments(self,tmp_path,default):
         default.add_comments(['this is a comment'])
         default.save(tmp_path/'with_comments',parallel=False)
-        new = VTK.load(tmp_path/'with_comments.vtr')
+        new = VTK.load(tmp_path/'with_comments.vti')
         assert new.get_comments() == ['this is a comment']
 
     @pytest.mark.xfail(int(vtk.vtkVersion.GetVTKVersion().split('.')[0])<8, reason='missing METADATA')

python/tests/test_util.py

@@ -59,9 +59,22 @@ class TestUtil:
                               ([1,1,0],'x',False,False,[0.5,0]),
                               ([1,1,1],'y',True, True, [0.3660254, 0,0.3660254]),
                              ])
-    def test_project_stereographic(self,point,direction,normalize,keepdims,answer):
+    def test_project_equal_angle(self,point,direction,normalize,keepdims,answer):
-        assert np.allclose(util.project_stereographic(np.array(point),direction=direction,
+        assert np.allclose(util.project_equal_angle(np.array(point),direction=direction,
-                                                      normalize=normalize,keepdims=keepdims),answer)
+                                                    normalize=normalize,keepdims=keepdims),answer)
+
+    @pytest.mark.parametrize('point,direction,normalize,keepdims,answer',
+                             [
+                              ([1,0,0],'z',False,True, [1,0,0]),
+                              ([1,0,0],'z',True, False,[1,0]),
+                              ([0,1,1],'z',False,True, [0,0.70710678,0]),
+                              ([0,1,1],'y',True, False,[0.5411961,0]),
+                              ([1,1,0],'x',False,False,[0.70710678,0]),
+                              ([1,1,1],'y',True, True, [0.45970084,0,0.45970084]),
+                             ])
+    def test_project_equal_area(self,point,direction,normalize,keepdims,answer):
+        assert np.allclose(util.project_equal_area(np.array(point),direction=direction,
+                                                   normalize=normalize,keepdims=keepdims),answer)
 
     @pytest.mark.parametrize('fro,to,mode,answer',
                              [

src/LAPACK_interface.f90

@@ -6,7 +6,7 @@
 module LAPACK_interface
   interface
 
-    subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
+    pure subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
       use prec
       character,   intent(in)                             :: jobvl,jobvr
       integer,     intent(in)                             :: n,lda,ldvl,ldvr,lwork
@@ -18,16 +18,16 @@ module LAPACK_interface
       integer,     intent(out)                            :: info
     end subroutine dgeev
 
-    subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
+    pure subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
       use prec
       integer,     intent(in)                             :: n,nrhs,lda,ldb
       real(pReal), intent(inout), dimension(lda,n)        :: a
       integer,     intent(out),   dimension(n)            :: ipiv
-      real(pReal), intent(out),   dimension(ldb,nrhs)     :: b
+      real(pReal), intent(inout), dimension(ldb,nrhs)     :: b
       integer,     intent(out)                            :: info
     end subroutine dgesv
 
-    subroutine dgetrf(m,n,a,lda,ipiv,info)
+    pure subroutine dgetrf(m,n,a,lda,ipiv,info)
       use prec
       integer,     intent(in)                             :: m,n,lda
       real(pReal), intent(inout), dimension(lda,n)        :: a
@@ -35,16 +35,16 @@ module LAPACK_interface
       integer,     intent(out)                            :: info
     end subroutine dgetrf
 
-    subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
+    pure subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
       use prec
       integer,     intent(in)                             :: n,lda,lwork
       real(pReal), intent(inout), dimension(lda,n)        :: a
-      integer,     intent(out),   dimension(n)            :: ipiv
+      integer,     intent(in),    dimension(n)            :: ipiv
       real(pReal), intent(out),   dimension(max(1,lwork)) :: work
       integer,     intent(out)                            :: info
     end subroutine dgetri
 
-    subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
+    pure subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
       use prec
       character,   intent(in)                             :: jobz,uplo
       integer,     intent(in)                             :: n,lda,lwork

src/constants.f90

@@ -9,7 +9,8 @@ module constants
   public
 
   real(pReal), parameter :: &
-    T_ROOM = 300.0_pReal, &                                                                         !< Room temperature in K
+    T_ROOM = 300.0_pReal, &                                                                         !< Room temperature in K. ToDo: IUPAC: 298.15
-    K_B = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
+    K_B = 1.38e-23_pReal, &                                                                         !< Boltzmann constant in J/Kelvin
+    N_A = 6.02214076e23_pReal                                                                      !< Avogadro constant in 1/mol
 
 end module constants

src/grid/DAMASK_grid.f90

@@ -43,6 +43,15 @@ program DAMASK_grid
     logical ::                  estimate_rate                                                       !< follow trajectory of former loadcase
   end type tLoadCase
 
+!--------------------------------------------------------------------------------------------------
+! field labels information
+  enum, bind(c); enumerator :: &
+    FIELD_UNDEFINED_ID, &
+    FIELD_MECH_ID, &
+    FIELD_THERMAL_ID, &
+    FIELD_DAMAGE_ID
+  end enum
+
   integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
 
 !--------------------------------------------------------------------------------------------------

src/grid/grid_damage_spectral.f90

@@ -167,7 +167,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
   integer :: i, j, k, ce
   type(tSolutionState) :: solution
   PetscInt  :: devNull
-  PetscReal :: phi_min, phi_max, stagNorm, solnNorm
+  PetscReal :: phi_min, phi_max, stagNorm
 
   PetscErrorCode :: ierr
   SNESConvergedReason :: reason
@@ -189,9 +189,8 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
     solution%iterationsNeeded = totalIter
   end if
   stagNorm = maxval(abs(phi_current - phi_stagInc))
-  solnNorm = maxval(abs(phi_current))
   call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr)
-  solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*solnNorm)
+  solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*maxval(phi_current))
   call MPI_Allreduce(MPI_IN_PLACE,solution%stagConverged,1,MPI_LOGICAL,MPI_LAND,MPI_COMM_WORLD,ierr)
   phi_stagInc = phi_current
 

src/grid/grid_thermal_spectral.f90

@@ -162,7 +162,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
   integer :: i, j, k, ce
   type(tSolutionState) :: solution
   PetscInt  :: devNull
-  PetscReal :: T_min, T_max, stagNorm, solnNorm
+  PetscReal :: T_min, T_max, stagNorm
 
   PetscErrorCode :: ierr
   SNESConvergedReason :: reason
@@ -184,9 +184,8 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
     solution%iterationsNeeded = totalIter
   end if
   stagNorm = maxval(abs(T_current - T_stagInc))
-  solnNorm = maxval(abs(T_current))
   call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr)
-  solution%stagConverged = stagNorm < max(num%eps_thermal_atol, num%eps_thermal_rtol*solnNorm)
+  solution%stagConverged = stagNorm < max(num%eps_thermal_atol, num%eps_thermal_rtol*maxval(T_current))
   call MPI_Allreduce(MPI_IN_PLACE,solution%stagConverged,1,MPI_LOGICAL,MPI_LAND,MPI_COMM_WORLD,ierr)
   T_stagInc = T_current
 

src/grid/spectral_utilities.f90

@@ -28,15 +28,6 @@ module spectral_utilities
 
   include 'fftw3-mpi.f03'
 
-!--------------------------------------------------------------------------------------------------
-! field labels information
-  enum, bind(c); enumerator :: &
-    FIELD_UNDEFINED_ID, &
-    FIELD_MECH_ID, &
-    FIELD_THERMAL_ID, &
-    FIELD_DAMAGE_ID
-  end enum
-
 !--------------------------------------------------------------------------------------------------
 ! grid related information information
   real(pReal), protected,  public                :: wgt                                             !< weighting factor 1/Nelems
@@ -139,11 +130,7 @@ module spectral_utilities
     utilities_calculateRate, &
     utilities_forwardField, &
     utilities_updateCoords, &
-    utilities_saveReferenceStiffness, &
+    utilities_saveReferenceStiffness
-    FIELD_UNDEFINED_ID, &
-    FIELD_MECH_ID, &
-    FIELD_THERMAL_ID, &
-    FIELD_DAMAGE_ID
 
 contains
 
@@ -388,21 +375,24 @@ subroutine utilities_updateGamma(C)
     gamma_hat =  cmplx(0.0_pReal,0.0_pReal,pReal)                                                   ! for the singular point and any non invertible A
     do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red
       if (any([i,j,k] /= 1)) then                                                                   ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
-        forall(l = 1:3, m = 1:3) &
+        do concurrent (l = 1:3, m = 1:3)
           xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset)
-        forall(l = 1:3, m = 1:3) &
+        end do
+        do concurrent(l = 1:3, m = 1:3)
           temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+        end do
         A(1:3,1:3) = temp33_complex%re; A(4:6,4:6) =  temp33_complex%re
         A(1:3,4:6) = temp33_complex%im; A(4:6,1:3) = -temp33_complex%im
         if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
           call math_invert(A_inv, err, A)
           temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
-          forall(l=1:3, m=1:3, n=1:3, o=1:3) &
+          do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
             gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* &
                                                    conjg(-xi1st(o,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset)
-        endif
+          end do
-      endif
+        end if
-    enddo; enddo; enddo
+      end if
+    end do; end do; end do
   endif
 
 end subroutine utilities_updateGamma
@@ -505,32 +495,37 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
   memoryEfficient: if (num%memory_efficient) then
     do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
       if (any([i,j,k+grid3Offset] /= 1)) then                                                       ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
-        forall(l = 1:3, m = 1:3) &
+        do concurrent(l = 1:3, m = 1:3)
           xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
-        forall(l = 1:3, m = 1:3) &
+        end do
+        do concurrent(l = 1:3, m = 1:3)
           temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+        end do
         A(1:3,1:3) = temp33_complex%re; A(4:6,4:6) =  temp33_complex%re
         A(1:3,4:6) = temp33_complex%im; A(4:6,1:3) = -temp33_complex%im
         if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
           call math_invert(A_inv, err, A)
           temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
-          forall(l=1:3, m=1:3, n=1:3, o=1:3) &
+          do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
             gamma_hat(l,m,n,o,1,1,1) =  temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k)
+          end do
         else
           gamma_hat(1:3,1:3,1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal)
-        endif
+        end if
-        forall(l = 1:3, m = 1:3) &
+        do concurrent(l = 1:3, m = 1:3)
           temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
+        end do
         tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
-      endif
+      end if
-    enddo; enddo; enddo
+    end do; end do; end do
   else memoryEfficient
     do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid1Red
-      forall(l = 1:3, m = 1:3) &
+      do concurrent(l = 1:3, m = 1:3)
         temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k))
+      end do
       tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
-    enddo; enddo; enddo
+    end do; end do; end do
-  endif memoryEfficient
+  end if memoryEfficient
 
   if (grid3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
 

src/lattice.f90

@@ -2070,7 +2070,7 @@ end function getlabels
 !> @brief Equivalent Poisson's ratio (ν)
 !> @details https://doi.org/10.1143/JPSJ.20.635
 !--------------------------------------------------------------------------------------------------
-function lattice_equivalent_nu(C,assumption) result(nu)
+pure function lattice_equivalent_nu(C,assumption) result(nu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
   character(len=5),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
@@ -2103,7 +2103,7 @@ end function lattice_equivalent_nu
 !> @brief Equivalent shear modulus (μ)
 !> @details https://doi.org/10.1143/JPSJ.20.635
 !--------------------------------------------------------------------------------------------------
-function lattice_equivalent_mu(C,assumption) result(mu)
+pure function lattice_equivalent_mu(C,assumption) result(mu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
   character(len=5),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)

src/math.f90

@@ -512,7 +512,7 @@ end subroutine math_invert33
 !--------------------------------------------------------------------------------------------------
 !> @brief Inversion of symmetriced 3x3x3x3 matrix
 !--------------------------------------------------------------------------------------------------
-function math_invSym3333(A)
+pure function math_invSym3333(A)
 
   real(pReal),dimension(3,3,3,3)            :: math_invSym3333
 
@@ -538,7 +538,7 @@ end function math_invSym3333
 !--------------------------------------------------------------------------------------------------
 !> @brief invert quadratic matrix of arbitrary dimension
 !--------------------------------------------------------------------------------------------------
-subroutine math_invert(InvA, error, A)
+pure subroutine math_invert(InvA, error, A)
 
   real(pReal), dimension(:,:),                 intent(in)  :: A
   real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA
@@ -895,7 +895,7 @@ pure function math_33toVoigt6_stress(sigma) result(sigma_tilde)
 
 
   sigma_tilde = [sigma(1,1), sigma(2,2), sigma(3,3), &
-                 sigma(3,2), sigma(3,1), sigma(1,2)] 
+                 sigma(3,2), sigma(3,1), sigma(1,2)]
 
 end function math_33toVoigt6_stress
 
@@ -910,7 +910,7 @@ pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde)
 
 
   epsilon_tilde = [          epsilon(1,1),           epsilon(2,2),           epsilon(3,3), &
-                   2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)] 
+                   2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)]
 
 end function math_33toVoigt6_strain
 
@@ -961,45 +961,42 @@ pure function math_3333toVoigt66(m3333)
 end function math_3333toVoigt66
 
 
-
 !--------------------------------------------------------------------------------------------------
-!> @brief draw a random sample from Gauss variable
+!> @brief Draw a sample from a normal distribution.
+!> @details https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
+!>          https://masuday.github.io/fortran_tutorial/random.html
 !--------------------------------------------------------------------------------------------------
-real(pReal) function math_sampleGaussVar(mu, sigma, width)
+impure elemental subroutine math_normal(x,mu,sigma)
 
-  real(pReal), intent(in) ::            mu, &                                                       !< mean
+  real(pReal), intent(out) :: x
-                                        sigma                                                       !< standard deviation
+  real(pReal), intent(in), optional :: mu, sigma
-  real(pReal), intent(in), optional ::  width                                                       !< cut off as multiples of standard deviation
 
-  real(pReal), dimension(2) ::          rnd                                                         ! random numbers
+  real(pReal) :: sigma_, mu_
-  real(pReal) ::                        scatter, &                                                  ! normalized scatter around mean
+  real(pReal), dimension(2) :: rnd
-                                        width_
 
-  if (abs(sigma) < tol_math_check) then
+
-    math_sampleGaussVar = mu
+  if (present(mu)) then
+    mu_ = mu
   else
-    if (present(width)) then
+    mu_ = 0.0_pReal
-      width_ = width
+  end if
-    else
-      width_ = 3.0_pReal                                                                            ! use +-3*sigma as default scatter
-    endif
 
-    do
+  if (present(sigma)) then
-      call random_number(rnd)
+    sigma_ = sigma
-      scatter = width_ * (2.0_pReal * rnd(1) - 1.0_pReal)
+  else
-      if (rnd(2) <= exp(-0.5_pReal * scatter**2)) exit                                              ! test if scattered value is drawn
+    sigma_ = 1.0_pReal
-    enddo
+  end if
 
-    math_sampleGaussVar = scatter * sigma
+  call random_number(rnd)
-  endif
+  x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(2.0_pReal*PI*(1.0_pReal - rnd(2)))
 
-end function math_sampleGaussVar
+end subroutine math_normal
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief eigenvalues and eigenvectors of symmetric matrix
 !--------------------------------------------------------------------------------------------------
-subroutine math_eigh(w,v,error,m)
+pure subroutine math_eigh(w,v,error,m)
 
   real(pReal), dimension(:,:),                  intent(in)  :: m                                    !< quadratic matrix to compute eigenvectors and values of
   real(pReal), dimension(size(m,1)),            intent(out) :: w                                    !< eigenvalues
@@ -1024,7 +1021,7 @@ end subroutine math_eigh
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
 !> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
 !--------------------------------------------------------------------------------------------------
-subroutine math_eigh33(w,v,m)
+pure subroutine math_eigh33(w,v,m)
 
   real(pReal), dimension(3,3),intent(in)  :: m                                                      !< 3x3 matrix to compute eigenvectors and values of
   real(pReal), dimension(3),  intent(out) :: w                                                      !< eigenvalues
@@ -1117,7 +1114,7 @@ end function math_rotationalPart
 !> @brief Eigenvalues of symmetric matrix
 ! will return NaN on error
 !--------------------------------------------------------------------------------------------------
-function math_eigvalsh(m)
+pure function math_eigvalsh(m)
 
   real(pReal), dimension(:,:),                  intent(in)  :: m                                    !< symmetric matrix to compute eigenvalues of
   real(pReal), dimension(size(m,1))                         :: math_eigvalsh
@@ -1140,7 +1137,7 @@ end function math_eigvalsh
 !> but apparently more stable solution and has general LAPACK powered version for arbritrary sized
 !> matrices as fallback
 !--------------------------------------------------------------------------------------------------
-function math_eigvalsh33(m)
+pure function math_eigvalsh33(m)
 
   real(pReal), intent(in), dimension(3,3) :: m                                                      !< 3x3 symmetric matrix to compute eigenvalues of
   real(pReal), dimension(3) :: math_eigvalsh33,I
@@ -1434,6 +1431,28 @@ subroutine selfTest
   if (dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) &
     error stop 'math_LeviCivita'
 
+  normal_distribution: block
+    integer, parameter :: N = 1000000
+    real(pReal), dimension(:), allocatable :: r
+    real(pReal) :: mu, sigma
+
+    allocate(r(N))
+    call random_number(mu)
+    call random_number(sigma)
+
+    sigma = 1.0_pReal + sigma*5.0_pReal
+    mu = (mu-0.5_pReal)*10_pReal
+
+    call math_normal(r,mu,sigma)
+
+    if (abs(mu -sum(r)/real(N,pReal))>5.0e-2_pReal) &
+      error stop 'math_normal(mu)'
+
+    mu = sum(r)/real(N,pReal)
+    if (abs(sigma**2 -1.0_pReal/real(N-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) &
+      error stop 'math_normal(sigma)'
+  end block normal_distribution
+
 end subroutine selfTest
 
 end module math

src/phase.f90

@@ -155,7 +155,7 @@ module phase
       real(pReal), dimension(3,3) :: P
     end function phase_P
 
-    module function thermal_T(ph,en) result(T)
+    pure module function thermal_T(ph,en) result(T)
       integer, intent(in) :: ph,en
       real(pReal) :: T
     end function thermal_T

src/phase_mechanical.f90

@@ -5,17 +5,17 @@ submodule(phase) mechanical
 
 
   enum, bind(c); enumerator :: &
-    PLASTICITY_UNDEFINED_ID, &
+    PLASTIC_UNDEFINED_ID, &
-    PLASTICITY_NONE_ID, &
+    PLASTIC_NONE_ID, &
-    PLASTICITY_ISOTROPIC_ID, &
+    PLASTIC_ISOTROPIC_ID, &
-    PLASTICITY_PHENOPOWERLAW_ID, &
+    PLASTIC_PHENOPOWERLAW_ID, &
-    PLASTICITY_KINEHARDENING_ID, &
+    PLASTIC_KINEHARDENING_ID, &
-    PLASTICITY_DISLOTWIN_ID, &
+    PLASTIC_DISLOTWIN_ID, &
-    PLASTICITY_DISLOTUNGSTEN_ID, &
+    PLASTIC_DISLOTUNGSTEN_ID, &
-    PLASTICITY_NONLOCAL_ID, &
+    PLASTIC_NONLOCAL_ID, &
-    KINEMATICS_UNDEFINED_ID, &
+    EIGEN_UNDEFINED_ID, &
-    KINEMATICS_CLEAVAGE_OPENING_ID, &
+    EIGEN_CLEAVAGE_OPENING_ID, &
-    KINEMATICS_THERMAL_EXPANSION_ID
+    EIGEN_THERMAL_EXPANSION_ID
   end enum
 
   type(tTensorContainer), dimension(:), allocatable :: &
@@ -37,7 +37,7 @@ submodule(phase) mechanical
     phase_mechanical_S0
 
 
-  integer(kind(PLASTICITY_undefined_ID)), dimension(:),   allocatable :: &
+  integer(kind(PLASTIC_undefined_ID)), dimension(:),   allocatable :: &
     phase_plasticity                                                                                !< plasticity of each phase
 
   integer :: phase_plasticity_maxSizeDotState
@@ -168,17 +168,17 @@ submodule(phase) mechanical
       integer,     intent(in) :: ph,en
     end function plastic_dislotwin_homogenizedC
 
-    module function elastic_C66(ph,en) result(C66)
+    pure module function elastic_C66(ph,en) result(C66)
       real(pReal), dimension(6,6) :: C66
       integer,     intent(in) :: ph, en
     end function elastic_C66
 
-    module function elastic_mu(ph,en) result(mu)
+    pure module function elastic_mu(ph,en) result(mu)
       real(pReal) :: mu
       integer, intent(in) :: ph, en
     end function elastic_mu
 
-    module function elastic_nu(ph,en) result(nu)
+    pure module function elastic_nu(ph,en) result(nu)
       real(pReal) :: nu
       integer, intent(in) :: ph, en
     end function elastic_nu
@@ -291,7 +291,7 @@ module subroutine mechanical_init(phases)
   call elastic_init(phases)
 
   allocate(plasticState(phases%length))
-  allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID)
+  allocate(phase_plasticity(phases%length),source = PLASTIC_UNDEFINED_ID)
   call plastic_init()
   do ph = 1,phases%length
     plasticState(ph)%state0 = plasticState(ph)%state
@@ -340,22 +340,22 @@ module subroutine mechanical_results(group,ph)
 
   select case(phase_plasticity(ph))
 
-    case(PLASTICITY_ISOTROPIC_ID)
+    case(PLASTIC_ISOTROPIC_ID)
       call plastic_isotropic_results(ph,group//'mechanical/')
 
-    case(PLASTICITY_PHENOPOWERLAW_ID)
+    case(PLASTIC_PHENOPOWERLAW_ID)
       call plastic_phenopowerlaw_results(ph,group//'mechanical/')
 
-    case(PLASTICITY_KINEHARDENING_ID)
+    case(PLASTIC_KINEHARDENING_ID)
       call plastic_kinehardening_results(ph,group//'mechanical/')
 
-    case(PLASTICITY_DISLOTWIN_ID)
+    case(PLASTIC_DISLOTWIN_ID)
       call plastic_dislotwin_results(ph,group//'mechanical/')
 
-    case(PLASTICITY_DISLOTUNGSTEN_ID)
+    case(PLASTIC_DISLOTUNGSTEN_ID)
       call plastic_dislotungsten_results(ph,group//'mechanical/')
 
-    case(PLASTICITY_NONLOCAL_ID)
+    case(PLASTIC_NONLOCAL_ID)
       call plastic_nonlocal_results(ph,group//'mechanical/')
 
   end select

src/phase_mechanical_eigen.f90

@@ -3,9 +3,9 @@ submodule(phase:mechanical) eigen
   integer, dimension(:), allocatable :: &
     Nmodels
 
-  integer(kind(KINEMATICS_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
+  integer(kind(EIGEN_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
     model
-  integer(kind(KINEMATICS_UNDEFINED_ID)),  dimension(:), allocatable :: &
+  integer(kind(EIGEN_UNDEFINED_ID)),  dimension(:), allocatable :: &
     model_damage
 
   interface
@@ -57,15 +57,15 @@ module subroutine eigen_init(phases)
     Nmodels(ph) = kinematics%length
   end do
 
-  allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID)
+  allocate(model(maxval(Nmodels),phases%length), source = EIGEN_undefined_ID)
 
   if (maxval(Nmodels) /= 0) then
-    where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID
+    where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID
   endif
 
-  allocate(model_damage(phases%length),  source = KINEMATICS_UNDEFINED_ID)
+  allocate(model_damage(phases%length),  source = EIGEN_UNDEFINED_ID)
 
-  where(damage_anisobrittle_init())  model_damage = KINEMATICS_cleavage_opening_ID
+  where(damage_anisobrittle_init())  model_damage = EIGEN_cleavage_opening_ID
 
 
 end subroutine eigen_init
@@ -173,7 +173,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
 
 
   plasticType: select case (phase_plasticity(ph))
-    case (PLASTICITY_isotropic_ID) plasticType
+    case (PLASTIC_isotropic_ID) plasticType
       call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en)
       Li = Li + my_Li
       dLi_dS = dLi_dS + my_dLi_dS
@@ -182,7 +182,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
 
   KinematicsLoop: do k = 1, Nmodels(ph)
     kinematicsType: select case (model(k,ph))
-      case (KINEMATICS_thermal_expansion_ID) kinematicsType
+      case (EIGEN_thermal_expansion_ID) kinematicsType
         call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en)
         Li = Li + my_Li
         dLi_dS = dLi_dS + my_dLi_dS
@@ -191,7 +191,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
   end do KinematicsLoop
 
   select case (model_damage(ph))
-    case (KINEMATICS_cleavage_opening_ID)
+    case (EIGEN_cleavage_opening_ID)
       call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en)
       Li = Li + my_Li
       dLi_dS = dLi_dS + my_dLi_dS

src/phase_mechanical_elastic.f90

@@ -86,7 +86,7 @@ end subroutine elastic_init
 !--------------------------------------------------------------------------------------------------
 !> @brief return 6x6 elasticity tensor
 !--------------------------------------------------------------------------------------------------
-module function elastic_C66(ph,en) result(C66)
+pure module function elastic_C66(ph,en) result(C66)
 
   integer, intent(in) :: &
     ph, &
@@ -140,7 +140,7 @@ end function elastic_C66
 !--------------------------------------------------------------------------------------------------
 !> @brief return shear modulus
 !--------------------------------------------------------------------------------------------------
-module function elastic_mu(ph,en) result(mu)
+pure module function elastic_mu(ph,en) result(mu)
 
   integer, intent(in) :: &
     ph, &
@@ -157,7 +157,7 @@ end function elastic_mu
 !--------------------------------------------------------------------------------------------------
 !> @brief return Poisson ratio
 !--------------------------------------------------------------------------------------------------
-module function elastic_nu(ph,en) result(nu)
+pure module function elastic_nu(ph,en) result(nu)
 
   integer, intent(in) :: &
     ph, &
@@ -223,7 +223,7 @@ module function phase_homogenizedC66(ph,en) result(C)
 
 
   plasticType: select case (phase_plasticity(ph))
-    case (PLASTICITY_DISLOTWIN_ID) plasticType
+    case (PLASTIC_DISLOTWIN_ID) plasticType
      C = plastic_dislotwin_homogenizedC(ph,en)
     case default plasticType
      C = elastic_C66(ph,en)

src/phase_mechanical_plastic.f90

@@ -73,47 +73,37 @@ submodule(phase:mechanical) plastic
         en
     end subroutine kinehardening_LpAndItsTangent
 
-    module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
+    module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Lp
       real(pReal), dimension(3,3,3,3), intent(out) :: &
         dLp_dMp
-
       real(pReal), dimension(3,3),     intent(in) :: &
         Mp
-      real(pReal),                     intent(in) :: &
-        T
       integer,                         intent(in) :: &
         ph, &
         en
     end subroutine dislotwin_LpAndItsTangent
 
-    pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
+    pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Lp
       real(pReal), dimension(3,3,3,3), intent(out) :: &
         dLp_dMp
-
       real(pReal), dimension(3,3),     intent(in) :: &
         Mp
-      real(pReal),                     intent(in) :: &
-        T
       integer,                         intent(in) :: &
         ph, &
         en
     end subroutine dislotungsten_LpAndItsTangent
 
-    module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
+    module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
-                                                       Mp,Temperature,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Lp
       real(pReal), dimension(3,3,3,3), intent(out) :: &
         dLp_dMp
-
       real(pReal), dimension(3,3),     intent(in) :: &
         Mp                                                                                          !< Mandel stress
-      real(pReal),                     intent(in) :: &
-        Temperature
       integer,                         intent(in) :: &
         ph, &
         en
@@ -224,15 +214,15 @@ module subroutine plastic_init
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic init  -+>>>'
 
-  where(plastic_none_init())              phase_plasticity = PLASTICITY_NONE_ID
+  where(plastic_none_init())              phase_plasticity = PLASTIC_NONE_ID
-  where(plastic_isotropic_init())         phase_plasticity = PLASTICITY_ISOTROPIC_ID
+  where(plastic_isotropic_init())         phase_plasticity = PLASTIC_ISOTROPIC_ID
-  where(plastic_phenopowerlaw_init())     phase_plasticity = PLASTICITY_PHENOPOWERLAW_ID
+  where(plastic_phenopowerlaw_init())     phase_plasticity = PLASTIC_PHENOPOWERLAW_ID
-  where(plastic_kinehardening_init())     phase_plasticity = PLASTICITY_KINEHARDENING_ID
+  where(plastic_kinehardening_init())     phase_plasticity = PLASTIC_KINEHARDENING_ID
-  where(plastic_dislotwin_init())         phase_plasticity = PLASTICITY_DISLOTWIN_ID
+  where(plastic_dislotwin_init())         phase_plasticity = PLASTIC_DISLOTWIN_ID
-  where(plastic_dislotungsten_init())     phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID
+  where(plastic_dislotungsten_init())     phase_plasticity = PLASTIC_DISLOTUNGSTEN_ID
-  where(plastic_nonlocal_init())          phase_plasticity = PLASTICITY_NONLOCAL_ID
+  where(plastic_nonlocal_init())          phase_plasticity = PLASTIC_NONLOCAL_ID
 
-  if (any(phase_plasticity == PLASTICITY_undefined_ID)) call IO_error(201)
+  if (any(phase_plasticity == PLASTIC_undefined_ID)) call IO_error(201)
 
 end subroutine plastic_init
 
@@ -262,7 +252,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
     i, j
 
 
-  if (phase_plasticity(ph) == PLASTICITY_NONE_ID) then
+  if (phase_plasticity(ph) == PLASTIC_NONE_ID) then
     Lp      = 0.0_pReal
     dLp_dFi = 0.0_pReal
     dLp_dS  = 0.0_pReal
@@ -272,23 +262,23 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
 
     plasticType: select case (phase_plasticity(ph))
 
-      case (PLASTICITY_ISOTROPIC_ID) plasticType
+      case (PLASTIC_ISOTROPIC_ID) plasticType
         call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTICITY_PHENOPOWERLAW_ID) plasticType
+      case (PLASTIC_PHENOPOWERLAW_ID) plasticType
         call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTICITY_KINEHARDENING_ID) plasticType
+      case (PLASTIC_KINEHARDENING_ID) plasticType
         call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTICITY_NONLOCAL_ID) plasticType
+      case (PLASTIC_NONLOCAL_ID) plasticType
-        call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en)
+        call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTICITY_DISLOTWIN_ID) plasticType
+      case (PLASTIC_DISLOTWIN_ID) plasticType
-        call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en)
+        call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
+      case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
-        call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en)
+        call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
     end select plasticType
 
@@ -321,28 +311,28 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
   logical :: broken
 
 
-  if (phase_plasticity(ph) /= PLASTICITY_NONE_ID) then
+  if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then
     Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
                        phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en))
 
     plasticType: select case (phase_plasticity(ph))
 
-      case (PLASTICITY_ISOTROPIC_ID) plasticType
+      case (PLASTIC_ISOTROPIC_ID) plasticType
         call isotropic_dotState(Mp,ph,en)
 
-      case (PLASTICITY_PHENOPOWERLAW_ID) plasticType
+      case (PLASTIC_PHENOPOWERLAW_ID) plasticType
         call phenopowerlaw_dotState(Mp,ph,en)
 
-      case (PLASTICITY_KINEHARDENING_ID) plasticType
+      case (PLASTIC_KINEHARDENING_ID) plasticType
         call plastic_kinehardening_dotState(Mp,ph,en)
 
-      case (PLASTICITY_DISLOTWIN_ID) plasticType
+      case (PLASTIC_DISLOTWIN_ID) plasticType
         call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en)
 
-      case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
+      case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
         call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en)
 
-      case (PLASTICITY_NONLOCAL_ID) plasticType
+      case (PLASTIC_NONLOCAL_ID) plasticType
         call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el)
     end select plasticType
   end if
@@ -372,13 +362,13 @@ module subroutine plastic_dependentState(co, ip, el)
 
   plasticType: select case (phase_plasticity(ph))
 
-    case (PLASTICITY_DISLOTWIN_ID) plasticType
+    case (PLASTIC_DISLOTWIN_ID) plasticType
       call dislotwin_dependentState(thermal_T(ph,en),ph,en)
 
-    case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
+    case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
       call dislotungsten_dependentState(ph,en)
 
-    case (PLASTICITY_NONLOCAL_ID) plasticType
+    case (PLASTIC_NONLOCAL_ID) plasticType
       call nonlocal_dependentState(ph,en,ip,el)
 
   end select plasticType
@@ -406,7 +396,7 @@ module function plastic_deltaState(ph, en) result(broken)
   broken = .false.
 
   select case (phase_plasticity(ph))
-    case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
+    case (PLASTIC_NONLOCAL_ID,PLASTIC_KINEHARDENING_ID)
 
       Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
                          phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
@@ -414,10 +404,10 @@ module function plastic_deltaState(ph, en) result(broken)
 
       plasticType: select case (phase_plasticity(ph))
 
-        case (PLASTICITY_KINEHARDENING_ID) plasticType
+        case (PLASTIC_KINEHARDENING_ID) plasticType
           call plastic_kinehardening_deltaState(Mp,ph,en)
 
-        case (PLASTICITY_NONLOCAL_ID) plasticType
+        case (PLASTIC_NONLOCAL_ID) plasticType
           call plastic_nonlocal_deltaState(Mp,ph,en)
 
       end select plasticType

src/phase_mechanical_plastic_dislotungsten.f90

@@ -257,27 +257,27 @@ end function plastic_dislotungsten_init
 !> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
 pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
-                                                     Mp,T,ph,en)
+                                                     Mp,ph,en)
   real(pReal), dimension(3,3),     intent(out) :: &
     Lp                                                                                              !< plastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out) :: &
     dLp_dMp                                                                                         !< derivative of Lp with respect to the Mandel stress
-
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< Mandel stress
-  real(pReal),                 intent(in) :: &
-    T                                                                                               !< temperature
   integer,                     intent(in) :: &
     ph, &
     en
 
   integer :: &
     i,k,l,m,n
+  real(pReal) :: &
+    T                                                                                               !< temperature
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     dot_gamma_pos,dot_gamma_neg, &
     ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
 
 
+  T = thermal_T(ph,en)
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
 

src/phase_mechanical_plastic_dislotwin.f90

@@ -476,18 +476,18 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
     C66_tw, &
     C66_tr
   integer :: i
-  real(pReal) :: f_unrotated
+  real(pReal) :: f_matrix
 
 
   C = elastic_C66(ph,en)
 
   associate(prm => param(ph), stt => state(ph))
 
-    f_unrotated = 1.0_pReal &
+    f_matrix = 1.0_pReal &
-                - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
+             - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
-                - sum(stt%f_tr(1:prm%sum_N_tr,en))
+             - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-    homogenizedC = f_unrotated * C
+    homogenizedC = f_matrix * C
 
     twinActive: if (prm%sum_N_tw > 0) then
       C66_tw    = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph))
@@ -513,20 +513,20 @@ end function plastic_dislotwin_homogenizedC
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
-module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
+module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
   real(pReal), dimension(3,3),     intent(out) :: Lp
   real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp
   real(pReal), dimension(3,3),     intent(in)  :: Mp
   integer,                         intent(in)  :: ph,en
-  real(pReal),                     intent(in)  :: T
 
   integer :: i,k,l,m,n
   real(pReal) :: &
-     f_unrotated,StressRatio_p,&
+    f_matrix,StressRatio_p,&
-     E_kB_T, &
+    E_kB_T, &
-     ddot_gamma_dtau, &
+    ddot_gamma_dtau, &
-     tau
+    tau, &
+    T
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     dot_gamma_sl,ddot_gamma_dtau_sl
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
@@ -556,69 +556,71 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
          0, 1, 1  &
          ],pReal),[ 3,6])
 
-  associate(prm => param(ph), stt => state(ph))
-
-  f_unrotated = 1.0_pReal &
-              - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
-              - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
+  T = thermal_T(ph,en)
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
 
-  call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl)
+  associate(prm => param(ph), stt => state(ph))
-  slipContribution: do i = 1, prm%sum_N_sl
-    Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
-    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
-  end do slipContribution
 
-  call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
+    f_matrix = 1.0_pReal &
-  twinContibution: do i = 1, prm%sum_N_tw
+             - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
-    Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
+             - sum(stt%f_tr(1:prm%sum_N_tr,en))
-    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
-  end do twinContibution
 
-  call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
+    call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl)
-  transContibution: do i = 1, prm%sum_N_tr
+    slipContribution: do i = 1, prm%sum_N_sl
-    Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
+      Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
-    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+        dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
+                         + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
-  end do transContibution
+    end do slipContribution
 
-  Lp      = Lp      * f_unrotated
+    call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
-  dLp_dMp = dLp_dMp * f_unrotated
+    twinContibution: do i = 1, prm%sum_N_tw
+      Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
+      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+        dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+                         + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
+    end do twinContibution
 
-  shearBandingContribution: if (dNeq0(prm%v_sb)) then
+    call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
+    transContibution: do i = 1, prm%sum_N_tr
+      Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
+      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+        dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+                         + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
+    end do transContibution
 
-    E_kB_T = prm%E_sb/(K_B*T)
+    Lp      = Lp      * f_matrix
-    call math_eigh33(eigValues,eigVectors,Mp)                                                       ! is Mp symmetric by design?
+    dLp_dMp = dLp_dMp * f_matrix
 
-    do i = 1,6
+    shearBandingContribution: if (dNeq0(prm%v_sb)) then
-      P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
-                                    matmul(eigVectors,sb_mComposition(1:3,i)))
-      tau = math_tensordot(Mp,P_sb)
 
-      significantShearBandStress: if (abs(tau) > tol_math_check) then
+      E_kB_T = prm%E_sb/(K_B*T)
-        StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
+      call math_eigh33(eigValues,eigVectors,Mp)                                                       ! is Mp symmetric by design?
-        dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
-        ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb &
-                        * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
-                        * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
 
-        Lp = Lp + dot_gamma_sb * P_sb
+      do i = 1,6
-        forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+        P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
-          dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+                                      matmul(eigVectors,sb_mComposition(1:3,i)))
-                           + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n)
+        tau = math_tensordot(Mp,P_sb)
-      end if significantShearBandStress
-    end do
 
-  end if shearBandingContribution
+        significantShearBandStress: if (abs(tau) > tol_math_check) then
+          StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
+          dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
+          ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb &
+                          * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
+                          * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
 
-  end associate
+          Lp = Lp + dot_gamma_sb * P_sb
+          forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+            dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+                             + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n)
+        end if significantShearBandStress
+      end do
+
+    end if shearBandingContribution
+
+    end associate
 
 end subroutine dislotwin_LpAndItsTangent
 
@@ -638,7 +640,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
 
   integer :: i
   real(pReal) :: &
-    f_unrotated, &
+    f_matrix, &
     d_hat, &
     v_cl, &                                                                                         !< climb velocity
     tau, &
@@ -661,9 +663,9 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
     mu = elastic_mu(ph,en)
     nu = elastic_nu(ph,en)
 
-    f_unrotated = 1.0_pReal &
+    f_matrix = 1.0_pReal &
-                - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
+             - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
-                - sum(stt%f_tr(1:prm%sum_N_tr,en))
+             - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
     call kinetics_sl(Mp,T,ph,en,dot_gamma_sl)
     dot%gamma_sl(:,en) = abs(dot_gamma_sl)
@@ -709,10 +711,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
                       - dot_rho_dip_climb
 
     call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
-    dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char
+    dot%f_tw(:,en) = f_matrix*dot_gamma_tw/prm%gamma_char
 
     call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
-    dot%f_tr(:,en) = f_unrotated*dot_gamma_tr
+    dot%f_tr(:,en) = f_matrix*dot_gamma_tr
 
   end associate
 

src/phase_mechanical_plastic_kinehardening.f90

@@ -86,6 +86,8 @@ module function plastic_kinehardening_init() result(myPlasticity)
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:kinehardening init  -+>>>'
   print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
+  print'(/,1x,a)', 'J.A. Wollmershauser et al., International Journal of Fatigue 36:181–193, 2012'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijfatigue.2011.07.008'
 
   phases => config_material%get('phase')
   allocate(param(phases%length))

src/phase_mechanical_plastic_nonlocal.f90

@@ -741,7 +741,7 @@ end subroutine nonlocal_dependentState
 !> @brief calculates plastic velocity gradient and its tangent
 !--------------------------------------------------------------------------------------------------
 module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
-                                                   Mp,Temperature,ph,en)
+                                                   Mp,ph,en)
   real(pReal), dimension(3,3), intent(out) :: &
     Lp                                                                                              !< plastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out) :: &
@@ -749,9 +749,6 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
   integer, intent(in) :: &
     ph, &
     en
-  real(pReal), intent(in) :: &
-    Temperature                                                                                     !< temperature
-
   real(pReal), dimension(3,3), intent(in) :: &
     Mp
                                                                                                     !< derivative of Lp with respect to Mp
@@ -771,67 +768,72 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     tau, &                                                                                          !< resolved shear stress including backstress terms
     dot_gamma                                                                                       !< shear rate
+  real(pReal) :: &
+    Temperature                                                                                     !< temperature
+
+
+  Temperature = thermal_T(ph,en)
+  Lp = 0.0_pReal
+  dLp_dMp = 0.0_pReal
 
   associate(prm => param(ph),dst=>dependentState(ph),stt=>state(ph))
 
-  !*** shortcut to state variables
+    !*** shortcut to state variables
-  rho = getRho(ph,en)
+    rho = getRho(ph,en)
-  rhoSgl = rho(:,sgl)
+    rhoSgl = rho(:,sgl)
 
-  do s = 1,prm%sum_N_sl
+    do s = 1,prm%sum_N_sl
-    tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s))
+      tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s))
-    tauNS(s,1) = tau(s)
+      tauNS(s,1) = tau(s)
-    tauNS(s,2) = tau(s)
+      tauNS(s,2) = tau(s)
-    if (tau(s) > 0.0_pReal) then
+      if (tau(s) > 0.0_pReal) then
-      tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s))
+        tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s))
-      tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s))
+        tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s))
-    else
+      else
-      tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s))
+        tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s))
-      tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s))
+        tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s))
-    end if
+      end if
-  end do
-  tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
-  tau   = tau   + dst%tau_back(:,en)
-
-  ! edges
-  call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
-                tau, tauNS(:,1), dst%tau_pass(:,en),1,Temperature, ph)
-  v(:,2)         = v(:,1)
-  dv_dtau(:,2)   = dv_dtau(:,1)
-  dv_dtauNS(:,2) = dv_dtauNS(:,1)
-
-  !screws
-  if (prm%nonSchmidActive) then
-    do t = 3,4
-      call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
-                    tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph)
     end do
-  else
+    tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
-    v(:,3:4)         = spread(v(:,1),2,2)
+    tau   = tau   + dst%tau_back(:,en)
-    dv_dtau(:,3:4)   = spread(dv_dtau(:,1),2,2)
-    dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2)
-  end if
 
-  stt%v(:,en) = pack(v,.true.)
+    ! edges
+    call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
+                  tau, tauNS(:,1), dst%tau_pass(:,en),1,Temperature, ph)
+    v(:,2)         = v(:,1)
+    dv_dtau(:,2)   = dv_dtau(:,1)
+    dv_dtauNS(:,2) = dv_dtauNS(:,1)
 
-  !*** Bauschinger effect
+    !screws
-  forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
+    if (prm%nonSchmidActive) then
-    rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
+      do t = 3,4
+        call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
+                      tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph)
+      end do
+    else
+      v(:,3:4)         = spread(v(:,1),2,2)
+      dv_dtau(:,3:4)   = spread(dv_dtau(:,1),2,2)
+      dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2)
+    end if
 
-  dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
+    stt%v(:,en) = pack(v,.true.)
 
-  Lp = 0.0_pReal
+    !*** Bauschinger effect
-  dLp_dMp = 0.0_pReal
+    forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
-  do s = 1,prm%sum_N_sl
+      rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
-    Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s)
+
-    forall (i=1:3,j=1:3,k=1:3,l=1:3) &
+    dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
-      dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
+
-          + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) &
+    do s = 1,prm%sum_N_sl
-          * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
+      Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s)
-          + prm%P_sl(i,j,s) &
+      forall (i=1:3,j=1:3,k=1:3,l=1:3) &
-          * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
+        dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
-             - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4))  * prm%b_sl(s)
+            + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) &
-  end do
+            * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
+            + prm%P_sl(i,j,s) &
+            * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
+               - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4))  * prm%b_sl(s)
+    end do
 
   end associate
 
@@ -870,7 +872,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
     dUpper, &                                                                                       ! current maximum stable dipole distance for edges and screws
     dUpperOld, &                                                                                    ! old maximum stable dipole distance for edges and screws
     deltaDUpper                                                                                     ! change in maximum stable dipole distance for edges and screws
-   
+
+
   associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph))
 
   mu = elastic_mu(ph,en)
@@ -1277,7 +1280,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
       !* The entering flux from my neighbor will be distributed on my slip systems according to the
       !* compatibility
       if (neighbor_n > 0) then
-      if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID .and. &
+      if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTIC_NONLOCAL_ID .and. &
           any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then
 
         forall (s = 1:ns, t = 1:4)
@@ -1323,7 +1326,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
       !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
       !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
       if (opposite_n > 0) then
-      if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTICITY_NONLOCAL_ID) then
+      if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then
 
         normal_me2neighbor_defConf = math_det33(Favg) &
                                    * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el))  ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
@@ -1394,78 +1397,79 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
     belowThreshold
   type(rotation) :: mis
 
+
   associate(prm => param(ph))
-  ns = prm%sum_N_sl
+    ns = prm%sum_N_sl
 
-  en = material_phaseMemberAt(1,i,e)
+    en = material_phaseMemberAt(1,i,e)
-  !*** start out fully compatible
+    !*** start out fully compatible
-  my_compatibility = 0.0_pReal
+    my_compatibility = 0.0_pReal
-  forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
+    forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
 
-  neighbors: do n = 1,nIPneighbors
+    neighbors: do n = 1,nIPneighbors
-    neighbor_e = IPneighborhood(1,n,i,e)
+      neighbor_e = IPneighborhood(1,n,i,e)
-    neighbor_i = IPneighborhood(2,n,i,e)
+      neighbor_i = IPneighborhood(2,n,i,e)
-    neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
+      neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
-    neighbor_phase = material_phaseAt(1,neighbor_e)
+      neighbor_phase = material_phaseAt(1,neighbor_e)
 
-    if (neighbor_e <= 0 .or. neighbor_i <= 0) then
+      if (neighbor_e <= 0 .or. neighbor_i <= 0) then
-      !* FREE SURFACE
+        !* FREE SURFACE
-      forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface)
+        forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface)
-    elseif (neighbor_phase /= ph) then
+      elseif (neighbor_phase /= ph) then
-      !* PHASE BOUNDARY
+        !* PHASE BOUNDARY
-      if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) &
+        if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) &
-        forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal
+          forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal
-    elseif (prm%chi_GB >= 0.0_pReal) then
+      elseif (prm%chi_GB >= 0.0_pReal) then
-      !* GRAIN BOUNDARY
+        !* GRAIN BOUNDARY
-      if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), &
+        if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), &
-                   phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. &
+                     phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. &
-          plasticState(neighbor_phase)%nonlocal) &
+            plasticState(neighbor_phase)%nonlocal) &
-        forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
+          forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
-    else
+      else
-      !* GRAIN BOUNDARY ?
+        !* GRAIN BOUNDARY ?
-      !* Compatibility defined by relative orientation of slip systems:
+        !* Compatibility defined by relative orientation of slip systems:
-      !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection.
+        !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection.
-      !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection.
+        !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection.
-      !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one),
+        !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one),
-      !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that
+        !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that
-      !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
+        !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
-      !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
+        !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
-      !* All values below the threshold are set to zero.
+        !* All values below the threshold are set to zero.
-      mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me))
+        mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me))
-      mySlipSystems: do s1 = 1,ns
+        mySlipSystems: do s1 = 1,ns
-        neighborSlipSystems: do s2 = 1,ns
+          neighborSlipSystems: do s2 = 1,ns
-          my_compatibility(1,s2,s1,n) =     math_inner(prm%slip_normal(1:3,s1), &
+            my_compatibility(1,s2,s1,n) =     math_inner(prm%slip_normal(1:3,s1), &
-                                                       mis%rotate(prm%slip_normal(1:3,s2))) &
+                                                         mis%rotate(prm%slip_normal(1:3,s2))) &
-                                      * abs(math_inner(prm%slip_direction(1:3,s1), &
+                                        * abs(math_inner(prm%slip_direction(1:3,s1), &
-                                                       mis%rotate(prm%slip_direction(1:3,s2))))
+                                                         mis%rotate(prm%slip_direction(1:3,s2))))
-          my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
+            my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
-                                                       mis%rotate(prm%slip_normal(1:3,s2)))) &
+                                                         mis%rotate(prm%slip_normal(1:3,s2)))) &
-                                      * abs(math_inner(prm%slip_direction(1:3,s1), &
+                                        * abs(math_inner(prm%slip_direction(1:3,s1), &
-                                                       mis%rotate(prm%slip_direction(1:3,s2))))
+                                                         mis%rotate(prm%slip_direction(1:3,s2))))
-        end do neighborSlipSystems
+          end do neighborSlipSystems
 
-        my_compatibilitySum = 0.0_pReal
+          my_compatibilitySum = 0.0_pReal
-        belowThreshold = .true.
+          belowThreshold = .true.
-        do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold))
+          do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold))
-          thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold)                       ! screws always positive
+            thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold)                     ! screws always positive
-          nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal)
+            nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal)
-          where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false.
+            where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false.
-          if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
+            if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
-            where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) &
+              where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) &
-              my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,&
+                my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,&
-                                                 my_compatibility(:,:,s1,n))
+                                                   my_compatibility(:,:,s1,n))
-          my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
+            my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
-        end do
+          end do
 
-        where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal
+          where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal
-        where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal
+          where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal
 
-      end do mySlipSystems
+        end do mySlipSystems
-    end if
+      end if
 
-  end do neighbors
+    end do neighbors
 
-  compatibility(:,:,:,:,i,e) = my_compatibility
+    compatibility(:,:,:,:,i,e) = my_compatibility
 
   end associate
 
@@ -1592,21 +1596,15 @@ subroutine stateInit(ini,phase,Nentries)
         stt%rhoSglMobile(s,e) = densityBinning
       end do
     else                                ! homogeneous distribution with noise
-      do e = 1, Nentries
+      do f = 1,size(ini%N_sl,1)
-        do f = 1,size(ini%N_sl,1)
+        from = 1 + sum(ini%N_sl(1:f-1))
-          from = 1 + sum(ini%N_sl(1:f-1))
+        upto = sum(ini%N_sl(1:f))
-          upto = sum(ini%N_sl(1:f))
+        call math_normal(stt%rho_sgl_mob_edg_pos(from:upto,:),ini%rho_u_ed_pos_0(f),ini%sigma_rho_u)
-          do s = from,upto
+        call math_normal(stt%rho_sgl_mob_edg_neg(from:upto,:),ini%rho_u_ed_neg_0(f),ini%sigma_rho_u)
-            noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), &
+        call math_normal(stt%rho_sgl_mob_scr_pos(from:upto,:),ini%rho_u_sc_pos_0(f),ini%sigma_rho_u)
-                     math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)]
+        call math_normal(stt%rho_sgl_mob_scr_neg(from:upto,:),ini%rho_u_sc_neg_0(f),ini%sigma_rho_u)
-            stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1)
+        stt%rho_dip_edg(from:upto,:) = ini%rho_d_ed_0(f)
-            stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1)
+        stt%rho_dip_scr(from:upto,:) = ini%rho_d_sc_0(f)
-            stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2)
-            stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2)
-          end do
-          stt%rho_dip_edg(from:upto,e)   = ini%rho_d_ed_0(f)
-          stt%rho_dip_scr(from:upto,e)   = ini%rho_d_sc_0(f)
-        end do
       end do
     end if
 
@@ -1652,53 +1650,55 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T,
     criticalStress_P, &                                                                             !< maximum obstacle strength
     criticalStress_S                                                                                !< maximum obstacle strength
 
-  associate(prm => param(ph))
+
   v = 0.0_pReal
   dv_dtau = 0.0_pReal
   dv_dtauNS = 0.0_pReal
 
-  do s = 1,prm%sum_N_sl
+  associate(prm => param(ph))
-    if (abs(tau(s)) > tauThreshold(s)) then
 
-      !* Peierls contribution
+    do s = 1,prm%sum_N_sl
-      tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s))
+      if (abs(tau(s)) > tauThreshold(s)) then
-      lambda_P = prm%b_sl(s)
-      activationVolume_P = prm%w *prm%b_sl(s)**3
-      criticalStress_P = prm%peierlsStress(s,c)
-      activationEnergy_P = criticalStress_P * activationVolume_P
-      tauRel_P = min(1.0_pReal, tauEff / criticalStress_P)
-      tPeierls = 1.0_pReal / prm%nu_a &
-               * exp(activationEnergy_P / (K_B * T) &
-                     * (1.0_pReal - tauRel_P**prm%p)**prm%q)
-      dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) &
-                             * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), &
-                             0.0_pReal, &
-                             tauEff < criticalStress_P)
 
-      ! Contribution from solid solution strengthening
+        !* Peierls contribution
-      tauEff = abs(tau(s)) - tauThreshold(s)
+        tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s))
-      lambda_S = prm%b_sl(s) / sqrt(prm%c_sol)
+        lambda_P = prm%b_sl(s)
-      activationVolume_S = prm%f_sol * prm%b_sl(s)**3 / sqrt(prm%c_sol)
+        activationVolume_P = prm%w *prm%b_sl(s)**3
-      criticalStress_S = prm%Q_sol / activationVolume_S
+        criticalStress_P = prm%peierlsStress(s,c)
-      tauRel_S = min(1.0_pReal, tauEff / criticalStress_S)
+        activationEnergy_P = criticalStress_P * activationVolume_P
-      tSolidSolution = 1.0_pReal /  prm%nu_a &
+        tauRel_P = min(1.0_pReal, tauEff / criticalStress_P)
-                     * exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
+        tPeierls = 1.0_pReal / prm%nu_a &
-      dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) &
+                 * exp(activationEnergy_P / (K_B * T) &
-                                   * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), &
+                       * (1.0_pReal - tauRel_P**prm%p)**prm%q)
-                                   0.0_pReal, &
+        dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) &
-                                   tauEff < criticalStress_S)
+                               * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), &
+                               0.0_pReal, &
+                               tauEff < criticalStress_P)
 
-      !* viscous glide velocity
+        ! Contribution from solid solution strengthening
-      tauEff = abs(tau(s)) - tauThreshold(s)
+        tauEff = abs(tau(s)) - tauThreshold(s)
+        lambda_S = prm%b_sl(s) / sqrt(prm%c_sol)
+        activationVolume_S = prm%f_sol * prm%b_sl(s)**3 / sqrt(prm%c_sol)
+        criticalStress_S = prm%Q_sol / activationVolume_S
+        tauRel_S = min(1.0_pReal, tauEff / criticalStress_S)
+        tSolidSolution = 1.0_pReal /  prm%nu_a &
+                       * exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
+        dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) &
+                                     * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), &
+                                     0.0_pReal, &
+                                     tauEff < criticalStress_S)
+
+        !* viscous glide velocity
+        tauEff = abs(tau(s)) - tauThreshold(s)
 
 
-      v(s) = sign(1.0_pReal,tau(s)) &
+        v(s) = sign(1.0_pReal,tau(s)) &
-           / (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff))
+             / (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff))
-      dv_dtau(s)   = v(s)**2 * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2))
+        dv_dtau(s)   = v(s)**2 * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2))
-      dv_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P
+        dv_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P
 
-    end if
+      end if
-  end do
+    end do
 
   end associate
 

src/phase_thermal.f90

@@ -271,7 +271,7 @@ end subroutine thermal_forward
 !----------------------------------------------------------------------------------------------
 !< @brief Get temperature (for use by non-thermal physics)
 !----------------------------------------------------------------------------------------------
-module function thermal_T(ph,en) result(T)
+pure module function thermal_T(ph,en) result(T)
 
   integer, intent(in) :: ph, en
   real(pReal) :: T