diff --git a/.github/workflows/Fortran.yml b/.github/workflows/Fortran.yml index 0d11918bc..5cc241e00 100644 --- a/.github/workflows/Fortran.yml +++ b/.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 - HOMEBREW_NO_AUTO_UPDATE: "ON" - HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: "ON" - HOMEBREW_NO_GITHUB_API: "ON" - HOMEBREW_NO_INSTALL_CLEANUP: "ON" + PETSC_VERSION: '3.16.2' + HOMEBREW_NO_ANALYTICS: 'ON' # Make Homebrew installation a little quicker + HOMEBREW_NO_AUTO_UPDATE: 'ON' + HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: '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 . - export PETSC_DIR=${PWD}/petsc-3.16.0 + tar -xf download/petsc-${PETSC_VERSION}.tar.gz -C . + 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 - key: petsc-3.16.0-${{ matrix.os }}-gcc${{ matrix.gcc_v }}-${{ hashFiles('**/petscversion.h') }} + path: petsc-${{ env.PETSC_VERSION }} + 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 . - export PETSC_DIR=${PWD}/petsc-3.16.0 + tar -xf download/petsc-${PETSC_VERSION}.tar.gz -C . + 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 - key: petsc-3.16.0-intel-${{ matrix.intel_v }}-${{ hashFiles('**/petscversion.h') }} + path: petsc-${{ env.PETSC_VERSION }} + 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 diff --git a/.github/workflows/Python.yml b/.github/workflows/Python.yml index c4fb15bdc..4c199b81c 100644 --- a/.github/workflows/Python.yml +++ b/.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 - python3-h5py python3-vtk7 python3-matplotlib python3-yaml -y + sudo apt-get remove mysql* && + 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 diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 756237b3f..6f18e3f1e 100644 --- a/.gitlab-ci.yml +++ b/.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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 56416d0fd..e8b774cfc 100644 --- a/CMakeLists.txt +++ b/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() diff --git a/PRIVATE b/PRIVATE index 02609955e..b898a8b55 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 02609955e53c0a5fb6cee5753419fb1ba1b9da2a +Subproject commit b898a8b5552bd9d1c555edc3d8134564dd32fe53 diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 5b551069e..3afffd2be 100644 --- a/cmake/Compiler-Intel.cmake +++ b/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 diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake new file mode 100644 index 000000000..326cfe319 --- /dev/null +++ b/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) diff --git a/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml b/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml index 3ad6779eb..f2066b373 100644 --- a/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml +++ b/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: - lattice: cI - mechanical: - elastic: {C_11: 417.4e+9, C_12: 242.4e+9, C_44: 211.1e+9, type: Hooke} - plastic: - N_sl: [12, 12] - a_sl: 2.0 - dot_gamma_0_sl: 0.001 - h_0_sl-sl: 563.0e+9 - 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] +N_sl: [12, 12] +a_sl: 2.0 +dot_gamma_0_sl: 0.001 +h_0_sl-sl: 563.0e+9 +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] diff --git a/examples/config/phase/AISI304.yaml b/examples/config/phase/AISI304.yaml new file mode 100644 index 000000000..0ca55635b --- /dev/null +++ b/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 diff --git a/examples/config/phase/Ag.yaml b/examples/config/phase/Ag.yaml new file mode 100644 index 000000000..ebcd21082 --- /dev/null +++ b/examples/config/phase/Ag.yaml @@ -0,0 +1,4 @@ +references: + - https://en.wikipedia.org/wiki/Silver +lattice: cF +rho: 10490.0 diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml new file mode 100644 index 000000000..6a3b339bb --- /dev/null +++ b/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 diff --git a/examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml b/examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml new file mode 100644 index 000000000..70fbd25c7 --- /dev/null +++ b/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 diff --git a/examples/config/phase/mechanical/elastic/Hooke_Ag.yaml b/examples/config/phase/mechanical/elastic/Hooke_Ag.yaml new file mode 100644 index 000000000..9485a4bb2 --- /dev/null +++ b/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 diff --git a/examples/config/phase/mechanical/elastic/Hooke_Al.yaml b/examples/config/phase/mechanical/elastic/Hooke_Al.yaml index 6007f7695..f4c266d54 100644 --- a/examples/config/phase/mechanical/elastic/Hooke_Al.yaml +++ b/examples/config/phase/mechanical/elastic/Hooke_Al.yaml @@ -1,8 +1,22 @@ type: Hooke references: - - J. Vallin et al., - Journal of Applied Physics 35(6):1825-1826, 1964, - https://doi.org/10.1063/1.1713749 -C_11: 107.3e+9 -C_12: 60.8e+9 -C_44: 28.3e+9 + - G.N. Kamm and G.A. Alers, + Journal of Applied Physics 35:327-330, 1964, + https://doi.org/10.1063/1.1713309 + - D. Gerlich and E.S. Fisher, + Journal of Physics and Chemistry of Solids 30:1197-1205, 1969 + 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 diff --git a/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml b/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml index e736829a6..94e761aaf 100644 --- a/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml +++ b/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml @@ -1,9 +1,19 @@ type: Hooke references: - - J.P. Hirth and J. Lothe, - Theory of Dislocations, 1982, - John Wiley & Sons, - page 837 -C_11: 242.e9 -C_12: 146.5e+9 -C_44: 112.e9 + - D.J. Dever, + Journal of Applied Physics 43(8):3293-3301, 1972, + https://doi.org/10.1063/1.1661710 + +T_ref: 300 + +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 diff --git a/examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml b/examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml new file mode 100644 index 000000000..7ba941066 --- /dev/null +++ b/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 diff --git a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml index 362797d0b..4e4ff1a9f 100644 --- a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml +++ b/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] diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml index 9b06f92d9..0ae04928d 100644 --- a/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml @@ -6,11 +6,11 @@ references: output: [xi_sl, xi_tw] -N_sl: [3, 3, 0, 6, 0, 6] # basal, prism, -, 1. pyr, -, 2. pyr -N_tw: [6, 0, 6] # tension, -, compression +N_sl: [3, 3, 6, 0, 6] # basal, prism, -, 1. pyr, -, 2. pyr +N_tw: [6, 0, 6] # tension, -, compression -xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6] -xi_inf_sl: [40.e+6, 135.e+6, 0., 150.e+6, 0., 150.e+6] +xi_0_sl: [10.e+6, 55.e+6, 60.e+6, 0., 60.e+6] +xi_inf_sl: [40.e+6, 135.e+6, 150.e+6, 0., 150.e+6] xi_0_tw: [40.e+6, 0., 60.e+6] a_sl: 2.25 @@ -23,16 +23,12 @@ f_sat_sl-tw: 10.0 h_0_sl-sl: 500.0e+6 h_0_tw-tw: 50.0e+6 h_0_tw-sl: 150.0e+6 -h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, - -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, - -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, - +1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, - +1.0, 1.0] # unused entries are indicated by -1.0 +h_sl-sl: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, + 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, +1.0, 1.0] # unused entries are indicated by -1.0 h_tw-tw: [+1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0] # unused entries are indicated by -1.0 -h_tw-sl: [+1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, - -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, - +1.0, -1.0, 1.0, -1.0] # unused entries are indicated by -1.0 -h_sl-tw: [+1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, - -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, - +1.0, -1.0, 1.0] # unused entries are indicated by -1.0 +h_tw-sl: [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, + 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, +1.0, -1.0, 1.0, -1.0] # unused entries are indicated by -1.0 +h_sl-tw: [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, + 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, +1.0, -1.0, 1.0] # unused entries are indicated by -1.0 diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml index 85fd2ee8e..1db6adc6b 100644 --- a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml @@ -8,18 +8,17 @@ references: https://doi.org/10.1016/j.actamat.2017.05.015 output: [gamma_sl] -N_sl: [3, 3, 0, 0, 12] # basal, 1. prism, -, -, 1. pyr +N_sl: [3, 3, 0, 12] # basal, 1. prism, -, 1. pyr n_sl: 20 a_sl: 2.0 dot_gamma_0_sl: 0.001 h_0_sl-sl: 200.e+6 # C. Zambaldi et al.: -xi_0_sl: [349.e+6, 150.e+6, 0.0, 0.0, 1107.e+6] -xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6] +xi_0_sl: [349.e+6, 150.e+6, 0.0, 1107.e+6] +xi_inf_sl: [568.e+6, 150.e+7, 0.0, 3420.e+6] # L. Wang et al. : -# xi_0_sl: [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6] +# xi_0_sl: [127.e+6, 96.e+6, 0.0, 240.e+6] h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, - -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, - +1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0] # unused entries are indicated by -1.0 + -1.0, -1.0, +1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0] # unused entries are indicated by -1.0 diff --git a/examples/config/phase/thermal/AISI304.yaml b/examples/config/phase/thermal/AISI304.yaml new file mode 100644 index 000000000..ef09fd850 --- /dev/null +++ b/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 diff --git a/examples/config/phase/thermal/Steel-0.5C.yaml b/examples/config/phase/thermal/steel-0.5C.yaml similarity index 100% rename from examples/config/phase/thermal/Steel-0.5C.yaml rename to examples/config/phase/thermal/steel-0.5C.yaml diff --git a/python/damask/VERSION b/python/damask/VERSION index 656df5ec5..a467f728c 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-244-gb57351045 +v3.0.0-alpha5-355-gc29428a60 diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 233be569a..c2721e0fa 100644 --- a/python/damask/_colormap.py +++ b/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], - high: Sequence[float], + def from_range(low: ArrayLike, + 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 \ - (min(bounds[:2]),max(bounds[:2])) + l,r = (field[mask].min(),field[mask].max()) if bounds is None else \ + 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 - hi,lo = hi+0.5*avg,lo-0.5*avg # extend range to have actual data centered within + if abs(delta) * 1e8 <= avg: # delta is similar to numerical noise + 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([ diff --git a/python/damask/_crystal.py b/python/damask/_crystal.py index fd87a90d4..689398cf6 100644 --- a/python/damask/_crystal.py +++ b/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. @@ -485,10 +487,6 @@ class Crystal(): [+2,-1,-1,+0, +0,+1,-1,+0], [-1,+2,-1,+0, -1,+0,+1,+0], [-1,-1,+2,+0, +1,-1,+0,+0]]), - np.array([ - [-1,+1,+0,+0, +1,+1,-2,+0], - [+0,-1,+1,+0, -2,+1,+1,+0], - [+1,+0,-1,+0, +1,-2,+1,+0]]), np.array([ [-1,+2,-1,+0, +1,+0,-1,+1], [-2,+1,+1,+0, +0,+1,-1,+1], @@ -555,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. diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index bb14fd38b..e727c54ae 100644 --- a/python/damask/_orientation.py +++ b/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 diff --git a/python/damask/_result.py b/python/damask/_result.py index 141b795c2..f47a7da80 100644 --- a/python/damask/_result.py +++ b/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_first = r_empty.view_more('increments',0) - >>> r_first_and_last = r.first.view_more('increments',-1) + >>> r_empty = damask.Result('my_file.hdf5').view(increments=False) + >>> r_first = r_empty.view_more(increments=0) + >>> 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 diff --git a/python/damask/util.py b/python/damask/util.py index a39a7b8c7..0581302db 100644 --- a/python/damask/util.py +++ b/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. - Default two-dimensional output uses right-handed frame spanned by + 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 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, - shift,axis=-1) - return np.roll(np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])),np.zeros_like(v_[...,2:3])]), + 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]/(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] diff --git a/python/tests/reference/Orientation/hP_slip.txt b/python/tests/reference/Orientation/hP_slip.txt index a67eea150..47896473b 100644 --- a/python/tests/reference/Orientation/hP_slip.txt +++ b/python/tests/reference/Orientation/hP_slip.txt @@ -5,9 +5,6 @@ 0.0 1.0 -5.914589856893347e-17 0.0 0.0 0.0 0.0 0.0 0.0 0.4330127018922192 0.24999999999999975 -3.102315069664884e-17 -0.7500000000000002 -0.4330127018922192 5.373367321746164e-17 0.0 0.0 0.0 -0.43301270189221935 0.25000000000000006 1.4502014121821253e-18 -0.7499999999999998 0.4330127018922193 2.5118225271075755e-18 0.0 0.0 0.0 --0.4330127018922194 -0.7499999999999999 6.059609998111558e-17 0.2500000000000001 0.4330127018922194 -3.498517463593857e-17 0.0 0.0 0.0 -2.563950248511418e-16 -5.693113199781536e-32 -9.614043519462407e-33 1.0 -2.220446049250313e-16 -3.7496997163046135e-17 0.0 0.0 0.0 -0.4330127018922194 -0.75 2.8122747872284606e-17 0.25000000000000006 -0.43301270189221935 1.6236676054415494e-17 0.0 0.0 0.0 -0.38254602783800284 -0.22086305214969287 -0.23426064283290896 0.6625891564490795 0.38254602783800284 0.40575133560034454 0.0 0.0 0.0 0.0 -0.8834522085987724 -0.4685212856658182 0.0 0.0 0.0 0.0 0.0 0.0 0.38254602783800307 -0.22086305214969315 -0.23426064283290912 0.6625891564490792 -0.38254602783800296 -0.40575133560034443 0.0 0.0 0.0 diff --git a/python/tests/test_Colormap.py b/python/tests/test_Colormap.py index ab9bcf92f..156907670 100644 --- a/python/tests/test_Colormap.py +++ b/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))) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 3f2555aec..4fea08eb3 100644 --- a/python/tests/test_Result.py +++ b/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',default.increments_in_range(0,np.iinfo(int).max)).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('times',True).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=True).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')) @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) diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index 26d1c4a53..68224ff33 100644 --- a/python/tests/test_VTK.py +++ b/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') diff --git a/python/tests/test_util.py b/python/tests/test_util.py index 13731eda2..ee345605f 100644 --- a/python/tests/test_util.py +++ b/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): - assert np.allclose(util.project_stereographic(np.array(point),direction=direction, - normalize=normalize,keepdims=keepdims),answer) + def test_project_equal_angle(self,point,direction,normalize,keepdims,answer): + assert np.allclose(util.project_equal_angle(np.array(point),direction=direction, + 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', [ diff --git a/src/LAPACK_interface.f90 b/src/LAPACK_interface.f90 index 7d3043ed0..e11f74875 100644 --- a/src/LAPACK_interface.f90 +++ b/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 diff --git a/src/constants.f90 b/src/constants.f90 index dd26ce78c..43bb60953 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -9,7 +9,8 @@ module constants public real(pReal), parameter :: & - T_ROOM = 300.0_pReal, & !< Room temperature in K - K_B = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin + T_ROOM = 300.0_pReal, & !< Room temperature in K. ToDo: IUPAC: 298.15 + K_B = 1.38e-23_pReal, & !< Boltzmann constant in J/Kelvin + N_A = 6.02214076e-23_pReal !< Avogadro constant in 1/mol end module constants diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index de324e0f5..b1da0c2a2 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/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(:) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 187a6f552..ef229368e 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/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 diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 868d3101e..67c7ba1c3 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/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 diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 29777db3b..f40025039 100644 --- a/src/grid/spectral_utilities.f90 +++ b/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, & - FIELD_UNDEFINED_ID, & - FIELD_MECH_ID, & - FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID + utilities_saveReferenceStiffness 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 - endif - enddo; enddo; enddo + end do + end if + 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 - forall(l = 1:3, m = 1:3) & + end if + 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 - enddo; enddo; enddo + end if + 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 - endif memoryEfficient + end do; end do; end do + end if memoryEfficient if (grid3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) diff --git a/src/lattice.f90 b/src/lattice.f90 index 935a8c161..3089aa30b 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -200,7 +200,7 @@ module lattice !-------------------------------------------------------------------------------------------------- ! hexagonal (hP) integer, dimension(*), parameter :: & - HEX_NSLIPSYSTEM = [3, 3, 3, 6, 12, 6] !< # of slip systems per family for hex + HEX_NSLIPSYSTEM = [3, 3, 6, 12, 6] !< # of slip systems per family for hex integer, dimension(*), parameter :: & HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex @@ -219,10 +219,6 @@ module lattice 2, -1, -1, 0, 0, 1, -1, 0, & -1, 2, -1, 0, -1, 0, 1, 0, & -1, -1, 2, 0, 1, -1, 0, 0, & - ! <-11.0>{11.0}/2. order prismatic compound systems (plane normal independent of c/a-ratio) - -1, 1, 0, 0, 1, 1, -2, 0, & - 0, -1, 1, 0, -2, 1, 1, 0, & - 1, 0, -1, 0, 1, -2, 1, 0, & ! <-1-1.0>{-11.1}/1. order pyramidal systems (direction independent of c/a-ratio) -1, 2, -1, 0, 1, 0, -1, 1, & -2, 1, 1, 0, 0, 1, -1, 1, & @@ -753,45 +749,41 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result( integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: & HEX_INTERACTIONSLIPSLIP = reshape( [& - ! basal prism 2. prism 1. pyr 1. pyr 2. pyr - 1, 2, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! -----> acting (forest) - 2, 1, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! | basal - 2, 2, 1, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! | - ! v - 6, 6, 6, 4, 5, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & ! reacting (primary) - 6, 6, 6, 5, 4, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & ! prism - 6, 6, 6, 5, 5, 4, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & + ! basal prism 1. pyr 1. pyr 2. pyr + 1, 2, 2, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13,13,13,13,13,13,13, 21,21,21,21,21,21, & ! -----> acting (forest) + 2, 1, 2, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13,13,13,13,13,13,13, 21,21,21,21,21,21, & ! | basal + 2, 2, 1, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13,13,13,13,13,13,13, 21,21,21,21,21,21, & ! | + ! v + 6, 6, 6, 4, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14,14,14,14,14,14,14, 22,22,22,22,22,22, & ! reacting (primary) + 6, 6, 6, 5, 4, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14,14,14,14,14,14,14, 22,22,22,22,22,22, & ! prism + 6, 6, 6, 5, 5, 4, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14,14,14,14,14,14,14, 22,22,22,22,22,22, & - 12,12,12, 11,11,11, 9,10,10, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, & - 12,12,12, 11,11,11, 10, 9,10, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, & ! 2. prism - 12,12,12, 11,11,11, 10,10, 9, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, & + 12,12,12, 11,11,11, 9,10,10,10,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & + 12,12,12, 11,11,11, 10, 9,10,10,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & + 12,12,12, 11,11,11, 10,10, 9,10,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & + 12,12,12, 11,11,11, 10,10,10, 9,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & ! 1. pyr + 12,12,12, 11,11,11, 10,10,10,10, 9,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & + 12,12,12, 11,11,11, 10,10,10,10,10, 9, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & - 20,20,20, 19,19,19, 18,18,18, 16,17,17,17,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & - 20,20,20, 19,19,19, 18,18,18, 17,16,17,17,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & - 20,20,20, 19,19,19, 18,18,18, 17,17,16,17,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & - 20,20,20, 19,19,19, 18,18,18, 17,17,17,16,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & ! 1. pyr - 20,20,20, 19,19,19, 18,18,18, 17,17,17,17,16,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & - 20,20,20, 19,19,19, 18,18,18, 17,17,17,17,17,16, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 16,17,17,17,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,16,17,17,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,16,17,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,16,17,17,17,17,17,17,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16,17,17,17,17,17,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,16,17,17,17,17,17, 24,24,24,24,24,24, & ! 1. pyr + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,16,17,17,17,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,16,17,17,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,17,16,17,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,17,17,16,17, 24,24,24,24,24,24, & + 20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,17,17,17,16, 24,24,24,24,24,24, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 25,26,26,26,26,26,26,26,26,26,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,25,26,26,26,26,26,26,26,26,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,25,26,26,26,26,26,26,26,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,25,26,26,26,26,26,26,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,25,26,26,26,26,26,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,25,26,26,26,26,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,25,26,26,26,26,26, 35,35,35,35,35,35, & ! 1. pyr - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,25,26,26,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,26,25,26,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,26,26,25,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,26,26,26,25,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,26,26,26,26,25, 35,35,35,35,35,35, & - - 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 36,37,37,37,37,37, & - 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,36,37,37,37,37, & - 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,36,37,37,37, & - 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,36,37,37, & ! 2. pyr - 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,36,37, & - 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,37,36 & + 30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 25,26,26,26,26,26, & + 30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,25,26,26,26,26, & + 30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,25,26,26,26, & + 30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,26,25,26,26, & ! 2. pyr + 30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,26,26,25,26, & + 30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,26,26,26,25 & ],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for hex (onion peel naming scheme) integer, dimension(BCT_NSLIP,BCT_NSLIP), parameter :: & @@ -1137,35 +1129,31 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & - 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & ! 2.prism + 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & + 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & + 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & ! 1. pyr + 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & - 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & ! 1. pyr + 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & + 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & + 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & + 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & ! 1. pyr + 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & + 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & + 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & + 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & ! 2. pyr 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & - 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & - 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & - 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & ! 1. pyr - 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & - 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & - 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & - 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & - 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & - - 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & - 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & - 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & - 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & ! 2. pyr - 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & - 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 & + 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20 & ],shape(HEX_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hex select case(lattice) @@ -1267,34 +1255,34 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) r integer, dimension(HEX_NSLIP,HEX_NTWIN), parameter :: & HEX_INTERACTIONTWINSLIP = reshape( [& - ! basal prism 2. prism 1. pyr 1. pyr 2. pyr - 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! ----> slip (acting) - 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! | - 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! | - 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! v <-10.1>{10.2} - 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! twin (reacting) - 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & + ! basal prism 1. pyr 1. pyr 2. pyr + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! ----> slip (acting) + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! | + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! | + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! v <-10.1>{10.2} + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! twin (reacting) + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & - 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & - 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & - 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & - 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & ! <11.6>{-1-1.1} - 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & - 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & + 2, 2, 2, 6, 6, 6, 10,10,10,10,10,10, 14,14,14,14,14,14,14,14,14,14,14,14, 18,18,18,18,18,18, & + 2, 2, 2, 6, 6, 6, 10,10,10,10,10,10, 14,14,14,14,14,14,14,14,14,14,14,14, 18,18,18,18,18,18, & + 2, 2, 2, 6, 6, 6, 10,10,10,10,10,10, 14,14,14,14,14,14,14,14,14,14,14,14, 18,18,18,18,18,18, & + 2, 2, 2, 6, 6, 6, 10,10,10,10,10,10, 14,14,14,14,14,14,14,14,14,14,14,14, 18,18,18,18,18,18, & ! <11.6>{-1-1.1} + 2, 2, 2, 6, 6, 6, 10,10,10,10,10,10, 14,14,14,14,14,14,14,14,14,14,14,14, 18,18,18,18,18,18, & + 2, 2, 2, 6, 6, 6, 10,10,10,10,10,10, 14,14,14,14,14,14,14,14,14,14,14,14, 18,18,18,18,18,18, & - 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & - 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & - 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & - 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & ! <10.-2>{10.1} - 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & - 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & + 3, 3, 3, 7, 7, 7, 11,11,11,11,11,11, 15,15,15,15,15,15,15,15,15,15,15,15, 19,19,19,19,19,19, & + 3, 3, 3, 7, 7, 7, 11,11,11,11,11,11, 15,15,15,15,15,15,15,15,15,15,15,15, 19,19,19,19,19,19, & + 3, 3, 3, 7, 7, 7, 11,11,11,11,11,11, 15,15,15,15,15,15,15,15,15,15,15,15, 19,19,19,19,19,19, & + 3, 3, 3, 7, 7, 7, 11,11,11,11,11,11, 15,15,15,15,15,15,15,15,15,15,15,15, 19,19,19,19,19,19, & ! <10.-2>{10.1} + 3, 3, 3, 7, 7, 7, 11,11,11,11,11,11, 15,15,15,15,15,15,15,15,15,15,15,15, 19,19,19,19,19,19, & + 3, 3, 3, 7, 7, 7, 11,11,11,11,11,11, 15,15,15,15,15,15,15,15,15,15,15,15, 19,19,19,19,19,19, & - 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & - 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & - 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & - 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & ! <11.-3>{11.2} - 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & - 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 & + 4, 4, 4, 8, 8, 8, 12,12,12,12,12,12, 16,16,16,16,16,16,16,16,16,16,16,16, 20,20,20,20,20,20, & + 4, 4, 4, 8, 8, 8, 12,12,12,12,12,12, 16,16,16,16,16,16,16,16,16,16,16,16, 20,20,20,20,20,20, & + 4, 4, 4, 8, 8, 8, 12,12,12,12,12,12, 16,16,16,16,16,16,16,16,16,16,16,16, 20,20,20,20,20,20, & + 4, 4, 4, 8, 8, 8, 12,12,12,12,12,12, 16,16,16,16,16,16,16,16,16,16,16,16, 20,20,20,20,20,20, & ! <11.-3>{11.2} + 4, 4, 4, 8, 8, 8, 12,12,12,12,12,12, 16,16,16,16,16,16,16,16,16,16,16,16, 20,20,20,20,20,20, & + 4, 4, 4, 8, 8, 8, 12,12,12,12,12,12, 16,16,16,16,16,16,16,16,16,16,16,16, 20,20,20,20,20,20 & ],shape(HEX_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hex select case(lattice) @@ -2082,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) @@ -2115,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) diff --git a/src/math.f90 b/src/math.f90 index 207a27e78..3ea4d6a08 100644 --- a/src/math.f90 +++ b/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 - sigma !< standard deviation - real(pReal), intent(in), optional :: width !< cut off as multiples of standard deviation + real(pReal), intent(out) :: x + real(pReal), intent(in), optional :: mu, sigma - real(pReal), dimension(2) :: rnd ! random numbers - real(pReal) :: scatter, & ! normalized scatter around mean - width_ + real(pReal) :: sigma_, mu_ + real(pReal), dimension(2) :: rnd - if (abs(sigma) < tol_math_check) then - math_sampleGaussVar = mu + + if (present(mu)) then + mu_ = mu else - if (present(width)) then - width_ = width - else - width_ = 3.0_pReal ! use +-3*sigma as default scatter - endif + mu_ = 0.0_pReal + end if - do - call random_number(rnd) - scatter = width_ * (2.0_pReal * rnd(1) - 1.0_pReal) - if (rnd(2) <= exp(-0.5_pReal * scatter**2)) exit ! test if scattered value is drawn - enddo + if (present(sigma)) then + sigma_ = sigma + else + sigma_ = 1.0_pReal + end if - math_sampleGaussVar = scatter * sigma - endif + call random_number(rnd) + 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 diff --git a/src/phase.f90 b/src/phase.f90 index 214b0f5fa..6035b4491 100644 --- a/src/phase.f90 +++ b/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 diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index db2faf914..1917d81c9 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -5,17 +5,17 @@ submodule(phase) mechanical enum, bind(c); enumerator :: & - PLASTICITY_UNDEFINED_ID, & - PLASTICITY_NONE_ID, & - PLASTICITY_ISOTROPIC_ID, & - PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_KINEHARDENING_ID, & - PLASTICITY_DISLOTWIN_ID, & - PLASTICITY_DISLOTUNGSTEN_ID, & - PLASTICITY_NONLOCAL_ID, & - KINEMATICS_UNDEFINED_ID, & - KINEMATICS_CLEAVAGE_OPENING_ID, & - KINEMATICS_THERMAL_EXPANSION_ID + PLASTIC_UNDEFINED_ID, & + PLASTIC_NONE_ID, & + PLASTIC_ISOTROPIC_ID, & + PLASTIC_PHENOPOWERLAW_ID, & + PLASTIC_KINEHARDENING_ID, & + PLASTIC_DISLOTWIN_ID, & + PLASTIC_DISLOTUNGSTEN_ID, & + PLASTIC_NONLOCAL_ID, & + EIGEN_UNDEFINED_ID, & + EIGEN_CLEAVAGE_OPENING_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 diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index c73665d17..f0b86319c 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/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 diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index ea113ddfb..24797a880 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/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) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index b2a9bdcc4..0a24b0cbf 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/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, & - Mp,Temperature,ph,en) + module subroutine nonlocal_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 !< 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_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID - where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTICITY_PHENOPOWERLAW_ID - where(plastic_kinehardening_init()) phase_plasticity = PLASTICITY_KINEHARDENING_ID - where(plastic_dislotwin_init()) phase_plasticity = PLASTICITY_DISLOTWIN_ID - where(plastic_dislotungsten_init()) phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID - where(plastic_nonlocal_init()) phase_plasticity = PLASTICITY_NONLOCAL_ID + where(plastic_none_init()) phase_plasticity = PLASTIC_NONE_ID + where(plastic_isotropic_init()) phase_plasticity = PLASTIC_ISOTROPIC_ID + where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTIC_PHENOPOWERLAW_ID + where(plastic_kinehardening_init()) phase_plasticity = PLASTIC_KINEHARDENING_ID + where(plastic_dislotwin_init()) phase_plasticity = PLASTIC_DISLOTWIN_ID + where(plastic_dislotungsten_init()) phase_plasticity = PLASTIC_DISLOTUNGSTEN_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 - call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + case (PLASTIC_NONLOCAL_ID) plasticType + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTICITY_DISLOTWIN_ID) plasticType - call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + case (PLASTIC_DISLOTWIN_ID) plasticType + call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + case (PLASTIC_DISLOTUNGSTEN_ID) plasticType + 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 diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index c2a584ea1..cd71a7fd7 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/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 diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index c5e043e3f..78681fa24 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/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 & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_matrix = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,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,& - E_kB_T, & - ddot_gamma_dtau, & - tau + f_matrix,StressRatio_p,& + E_kB_T, & + ddot_gamma_dtau, & + 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) - 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 + associate(prm => param(ph), stt => state(ph)) - call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) - 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 + f_matrix = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - 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 + call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) + 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 - Lp = Lp * f_unrotated - dLp_dMp = dLp_dMp * f_unrotated + call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) + 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) - call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? + Lp = Lp * f_matrix + dLp_dMp = dLp_dMp * f_matrix - do i = 1,6 - 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) + shearBandingContribution: if (dNeq0(prm%v_sb)) then - 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) + E_kB_T = prm%E_sb/(K_B*T) + call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? - 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 + do i = 1,6 + 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) - 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 & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_matrix = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,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 diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 6c0e1e0cc..03eb27f31 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/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)) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 6bd648eaa..7118055fe 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/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 - rho = getRho(ph,en) - rhoSgl = rho(:,sgl) + !*** shortcut to state variables + rho = getRho(ph,en) + rhoSgl = rho(:,sgl) - do s = 1,prm%sum_N_sl - tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) - tauNS(s,1) = tau(s) - tauNS(s,2) = tau(s) - if (tau(s) > 0.0_pReal) then - 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)) - else - 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)) - 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) + do s = 1,prm%sum_N_sl + tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + tauNS(s,1) = tau(s) + tauNS(s,2) = tau(s) + if (tau(s) > 0.0_pReal) then + 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)) + else + 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)) + end if 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 + tauNS = tauNS + spread(dst%tau_back(:,en),2,4) + tau = tau + dst%tau_back(:,en) - 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 - forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & - rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) + !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 + 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 - dLp_dMp = 0.0_pReal - do s = 1,prm%sum_N_sl - 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) & - dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & - + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) & - * 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 + !*** Bauschinger effect + forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & + rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) + + dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl + + do s = 1,prm%sum_N_sl + 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) & + dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & + + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) & + * 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) - !*** start out fully compatible - my_compatibility = 0.0_pReal - forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal + en = material_phaseMemberAt(1,i,e) + !*** start out fully compatible + my_compatibility = 0.0_pReal + forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal - neighbors: do n = 1,nIPneighbors - neighbor_e = IPneighborhood(1,n,i,e) - neighbor_i = IPneighborhood(2,n,i,e) - neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) - neighbor_phase = material_phaseAt(1,neighbor_e) + neighbors: do n = 1,nIPneighbors + neighbor_e = IPneighborhood(1,n,i,e) + neighbor_i = IPneighborhood(2,n,i,e) + neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) + neighbor_phase = material_phaseAt(1,neighbor_e) - if (neighbor_e <= 0 .or. neighbor_i <= 0) then - !* FREE SURFACE - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) - elseif (neighbor_phase /= ph) then - !* PHASE BOUNDARY - if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) & - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal - elseif (prm%chi_GB >= 0.0_pReal) then - !* GRAIN BOUNDARY - if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), & - phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & - plasticState(neighbor_phase)%nonlocal) & - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) - else - !* GRAIN BOUNDARY ? - !* 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. - !* 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), - !* 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. - !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. - !* All values below the threshold are set to zero. - mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) - mySlipSystems: do s1 = 1,ns - neighborSlipSystems: do s2 = 1,ns - my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & - mis%rotate(prm%slip_normal(1:3,s2))) & - * abs(math_inner(prm%slip_direction(1:3,s1), & - mis%rotate(prm%slip_direction(1:3,s2)))) - my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), & - mis%rotate(prm%slip_normal(1:3,s2)))) & - * abs(math_inner(prm%slip_direction(1:3,s1), & - mis%rotate(prm%slip_direction(1:3,s2)))) - end do neighborSlipSystems + if (neighbor_e <= 0 .or. neighbor_i <= 0) then + !* FREE SURFACE + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) + elseif (neighbor_phase /= ph) then + !* PHASE BOUNDARY + if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) & + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal + elseif (prm%chi_GB >= 0.0_pReal) then + !* GRAIN BOUNDARY + if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), & + phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & + plasticState(neighbor_phase)%nonlocal) & + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) + else + !* GRAIN BOUNDARY ? + !* 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. + !* 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), + !* 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. + !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. + !* All values below the threshold are set to zero. + mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) + mySlipSystems: do s1 = 1,ns + neighborSlipSystems: do s2 = 1,ns + my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & + mis%rotate(prm%slip_normal(1:3,s2))) & + * abs(math_inner(prm%slip_direction(1:3,s1), & + mis%rotate(prm%slip_direction(1:3,s2)))) + my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), & + mis%rotate(prm%slip_normal(1:3,s2)))) & + * abs(math_inner(prm%slip_direction(1:3,s1), & + mis%rotate(prm%slip_direction(1:3,s2)))) + end do neighborSlipSystems - my_compatibilitySum = 0.0_pReal - belowThreshold = .true. - do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold)) - thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold) ! screws always positive - nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal) - where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false. - if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) & - where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) & - my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,& - my_compatibility(:,:,s1,n)) - my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue - end do + my_compatibilitySum = 0.0_pReal + belowThreshold = .true. + do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold)) + thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold) ! screws always positive + nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal) + where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false. + if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) & + where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) & + my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,& + my_compatibility(:,:,s1,n)) + my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue + end do - where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal - where(belowThreshold) my_compatibility(2,:,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 - end do mySlipSystems - end if + end do mySlipSystems + 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) - from = 1 + sum(ini%N_sl(1:f-1)) - upto = sum(ini%N_sl(1:f)) - do s = from,upto - noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), & - math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)] - stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1) - stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) - 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 + do f = 1,size(ini%N_sl,1) + from = 1 + sum(ini%N_sl(1:f-1)) + 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) + call math_normal(stt%rho_sgl_mob_edg_neg(from:upto,:),ini%rho_u_ed_neg_0(f),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) + call math_normal(stt%rho_sgl_mob_scr_neg(from:upto,:),ini%rho_u_sc_neg_0(f),ini%sigma_rho_u) + stt%rho_dip_edg(from:upto,:) = ini%rho_d_ed_0(f) + stt%rho_dip_scr(from:upto,:) = ini%rho_d_sc_0(f) 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 - if (abs(tau(s)) > tauThreshold(s)) then + associate(prm => param(ph)) - !* Peierls contribution - tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) - 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) + do s = 1,prm%sum_N_sl + if (abs(tau(s)) > tauThreshold(s)) then - ! Contribution from solid solution strengthening - 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) + !* Peierls contribution + tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) + 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) - !* viscous glide velocity - tauEff = abs(tau(s)) - tauThreshold(s) + ! Contribution from solid solution strengthening + 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)) & - / (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_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P + v(s) = sign(1.0_pReal,tau(s)) & + / (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_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P - end if - end do + end if + end do end associate diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 4ea21d039..a3e4dd628 100644 --- a/src/phase_thermal.f90 +++ b/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