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

Merge branch 'development' into 23_BasticPETSc_2_PETSc

This commit is contained in:
Jaeyong Jung 2018-05-24 15:49:06 +02:00
parent 6800e779cb 35ba2604d4
commit c79a0d39c0
38 changed files with 1774 additions and 957 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
.gitlab-ci.yml
View File

@ -51,24 +51,24 @@ variables:
# ++++++++++++ Compiler ++++++++++++++++++++++++++++++++++++++++++++++
IntelCompiler16_0: "Compiler/Intel/16.0 Libraries/IMKL/2016"
IntelCompiler17_0: "Compiler/Intel/17.0 Libraries/IMKL/2017"
IntelCompiler18_1: "Compiler/Intel/18.1 Libraries/IMKL/2018"
GNUCompiler5_3: "Compiler/GNU/5.3"
# ------------ Defaults ----------------------------------------------
IntelCompiler: "$IntelCompiler17_0"
IntelCompiler: "$IntelCompiler18_1"
GNUCompiler: "$GNUCompiler5_3"
# ++++++++++++ MPI +++++++++++++++++++++++++++++++++++++++++++++++++++
MPICH3_2Intel17_0: "MPI/Intel/17.0/MPICH/3.2"
MPICH3_2GNU5_3: "MPI/GNU/5.3/MPICH/3.2"
MPICH3_2Intel18_1: "MPI/Intel/18.1/MPICH/3.2.1"
MPICH3_2GNU5_3: "MPI/GNU/5.3/MPICH/3.2.1"
# ------------ Defaults ----------------------------------------------
MPICH_GNU: "$MPICH3_2GNU5_3"
MPICH_Intel: "$MPICH3_2Intel17_0"
MPICH_Intel: "$MPICH3_2Intel18_1"
# ++++++++++++ PETSc +++++++++++++++++++++++++++++++++++++++++++++++++
PETSc3_7_6MPICH3_2Intel17_0: "Libraries/PETSc/3.7.6/Intel-17.0-MPICH-3.2"
PETSc3_7_5MPICH3_2Intel17_0: "Libraries/PETSc/3.7.5/Intel-17.0-MPICH-3.2"
PETSc3_6_4MPICH3_2Intel17_0: "Libraries/PETSc/3.6.4/Intel-17.0-MPICH-3.2"
PETSc3_7_5MPICH3_2GNU5_3: "Libraries/PETSc/3.7.5/GNU-5.3-MPICH-3.2"
PETSc3_9_1MPICH3_2Intel18_1: "Libraries/PETSc/3.9.1/Intel-18.1-MPICH-3.2.1"
PETSc3_9_1MPICH3_2GNU5_3: "Libraries/PETSc/3.9.1/GNU-5.3-MPICH-3.2.1"
# ------------ Defaults ----------------------------------------------
PETSc_MPICH_Intel: "$PETSc3_7_6MPICH3_2Intel17_0"
PETSc_MPICH_GNU: "$PETSc3_7_5MPICH3_2GNU5_3"
PETSc_MPICH_Intel: "$PETSc3_9_1MPICH3_2Intel18_1"
PETSc_MPICH_GNU: "$PETSc3_9_1MPICH3_2GNU5_3"
# ++++++++++++ FEM +++++++++++++++++++++++++++++++++++++++++++++++++++
Abaqus2016: "FEM/Abaqus/2016"
Abaqus2017: "FEM/Abaqus/2017"
@ -324,6 +324,13 @@ HybridIA:
- master
- release
TextureComponents:
stage: spectral
script: TextureComponents/test.py
except:
- master
- release
###################################################################################################
Marc_compileIfort2014:
stage: compileMarc2014

11
CMakeLists.txt
View File

@ -203,7 +203,9 @@ if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel")
# -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost"
endif ()
set (STANDARD_CHECK "-stand f08 -standard-semantics")
# -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 f08 -standard-semantics -assume nostd_mod_proc_name")
set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel")
# Link against shared Intel libraries instead of static ones
@ -215,13 +217,6 @@ if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel")
set (COMPILE_FLAGS "${COMPILE_FLAGS} -ftz")
# flush underflow to zero, automatically set if -O[1,2,3]
set (COMPILE_FLAGS "${COMPILE_FLAGS} -assume")
# assume ...
set (COMPILE_FLAGS "${COMPILE_FLAGS} byterecl")
# ... record length is given in bytes (also set by -standard-semantics)
set (COMPILE_FLAGS "${COMPILE_FLAGS},fpe_summary")
# ... print list of floating point exceptions occured during execution
set (COMPILE_FLAGS "${COMPILE_FLAGS} -diag-disable")
# disables warnings ...
set (COMPILE_FLAGS "${COMPILE_FLAGS} 5268")

1
DAMASK_env.csh
View File

@ -1 +0,0 @@
env/DAMASK.csh

1
DAMASK_env.sh
View File

@ -1 +0,0 @@
env/DAMASK.sh

1
DAMASK_env.zsh
View File

@ -1 +0,0 @@
env/DAMASK.zsh

2
DAMASK_prerequisites.sh
View File

@ -23,7 +23,7 @@ if which $1 &> /dev/null; then
$1 $2
echo -e '\n'
else
echo $ does not exist
echo $1 not found
fi
}

2
PRIVATE

@ -1 +1 @@
Subproject commit d3210425fcfd82673beffef5402049bd2954955c
Subproject commit bdcb4fdb26100a986c4cf224fe8bbd6f866934c9

2
VERSION
View File

@ -1 +1 @@
v2.0.1-1144-g5f28fa0
v2.0.1-1227-g7f0275e

6
installation/mods_Abaqus/abaqus_v6.env
View File

@ -26,14 +26,16 @@ fortCmd = "ifort"
# -ftz flush underflow to zero
# -diag-disable 5268 disable warnings about line length > 132 (only comments there anyway)
# -implicitnone assume no implicit types (e.g. i for integer)
# -assume byterecl count record length in bytes
# -standard-semantics sets standard (Fortran 2008) and some other conventions
# -assume nostd_mod_proc_name avoid problems with libraries compiled without that option
# -real-size 64 -DFLOAT=8 assume size of real to be 8 bytes, matches our definition of pReal
# -integer-size 32 -DINT=4 assume size of integer to be 4 bytes, matches our definition of pInt
compile_fortran = (fortCmd + " -c -fPIC -auto -shared-intel " +
"-I%I -free -O1 -fpp -openmp " +
"-ftz -diag-disable 5268 " +
"-implicitnone -assume byterecl -stand f08 -standard-semantics " +
"-implicitnone -standard-semantics " +
"-assume nostd_mod_proc_name " +
"-real-size 64 -integer-size 32 -DFLOAT=8 -DINT=4 " +
'-DDAMASKVERSION=\\\"%s\\\"'%DAMASKVERSION)

7
installation/mods_Abaqus/abaqus_v6_serial.env
View File

@ -19,21 +19,22 @@ fortCmd = "ifort"
# -free to use free-format FORTRAN 90 syntax
# -O <0-3> optimization level
# -fpp use FORTRAN preprocessor on source code
# -openmp build with openMP support
# -w90 -w95 suppress messages about use of non-standard Fortran (previous version of abaqus_v6.env only)
# -WB turn a compile-time bounds check into a warning (previous version of abaqus_v6.env only)
# -mP2OPT_hpo_vec_divbyzero=F inofficial compiler switch, proposed by abaqus but highly dubios (previous version of abaqus_v6.env only)
# -ftz flush underflow to zero
# -diag-disable 5268 disable warnings about line length > 132 (only comments there anyway)
# -implicitnone assume no implicit types (e.g. i for integer)
# -assume byterecl count record length in bytes
# -standard-semantics sets standard (Fortran 2008) and some other conventions
# -assume nostd_mod_proc_name avoid problems with libraries compiled without that option
# -real-size 64 -DFLOAT=8 assume size of real to be 8 bytes, matches our definition of pReal
# -integer-size 32 -DINT=4 assume size of integer to be 4 bytes, matches our definition of pInt
compile_fortran = (fortCmd + " -c -fPIC -auto -shared-intel " +
"-I%I -free -O1 -fpp " +
"-ftz -diag-disable 5268 " +
"-implicitnone -assume byterecl -stand f08 -standard-semantics " +
"-implicitnone -standard-semantics " +
"-assume nostd_mod_proc_name " +
"-real-size 64 -integer-size 32 -DFLOAT=8 -DINT=4 " +
'-DDAMASKVERSION=\\\"%s\\\"'%DAMASKVERSION)

14
installation/mods_MarcMentat/2014.2/Marc_tools/include_linux64
View File

@ -416,7 +416,7 @@ then
PROFILE=" $PROFILE -pg"
fi
FORT_OPT="-c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr -mp1 -WB"
FORT_OPT="-c -implicitnone -stand f08 -standard-semantics -assume nostd_std_mod_proc_name -safe_cray_ptr -mp1 -WB"
if test "$MTHREAD" = "OPENMP"
then
FORT_OPT=" $FORT_OPT -openmp"
@ -458,21 +458,21 @@ DAMASKVERSION="'"$DAMASKVERSION"'"
DFORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_std_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_std_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_std_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
@ -492,21 +492,21 @@ then
DFORTLOW="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_std_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2.2 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_std_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_std_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014.2 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"

14
installation/mods_MarcMentat/2014/Marc_tools/include_linux64
View File

@ -410,7 +410,7 @@ then
PROFILE="-prof-gen=srcpos"
fi
FORT_OPT="-c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr -mp1 -WB"
FORT_OPT="-c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr -mp1 -WB"
if test "$MTHREAD" = "OPENMP"
then
FORT_OPT=" $FORT_OPT -openmp"
@ -452,21 +452,21 @@ DAMASKVERSION="'"$DAMASKVERSION"'"
DFORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
@ -486,21 +486,21 @@ then
DFORTLOW="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2014 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"

14
installation/mods_MarcMentat/2015/Marc_tools/include_linux64
View File

@ -419,7 +419,7 @@ then
PROFILE=" $PROFILE -pg"
fi
FORT_OPT="-c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr -mp1 -WB"
FORT_OPT="-c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr -mp1 -WB"
if test "$MTHREAD" = "OPENMP"
then
FORT_OPT=" $FORT_OPT -openmp"
@ -454,21 +454,21 @@ fi
DFORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
@ -488,21 +488,21 @@ then
DFORTLOW="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2015 -DDAMASKVERSION=$DAMASKVERSION \
-openmp -openmp_report2 -openmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"

14
installation/mods_MarcMentat/2016/Marc_tools/include_linux64
View File

@ -449,7 +449,7 @@ then
PROFILE=" $PROFILE -pg"
fi
FORT_OPT="-c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr -mp1 -WB -fp-model source"
FORT_OPT="-c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr -mp1 -WB -fp-model source"
if test "$MTHREAD" = "OPENMP"
then
FORT_OPT=" $FORT_OPT -qopenmp"
@ -484,21 +484,21 @@ fi
DFORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
@ -518,21 +518,21 @@ then
DFORTLOW="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS"

14
installation/mods_MarcMentat/2017/Marc_tools/include_linux64
View File

@ -457,7 +457,7 @@ then
PROFILE=" $PROFILE -pg"
fi
FORT_OPT="-c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr -mp1 -WB -fp-model source"
FORT_OPT="-c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr -mp1 -WB -fp-model source"
if test "$MTHREAD" = "OPENMP"
then
FORT_OPT=" $FORT_OPT -qopenmp"
@ -494,21 +494,21 @@ fi
DFORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
@ -528,21 +528,21 @@ then
DFORTLOW="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTLOWMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRAN="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTRANMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGH="$FCOMP $FORT_OPT $PROFILE -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-enable sc3 -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
DFORTHIGHMP="$FCOMP -c -assume byterecl -stand f08 -standard-semantics -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
-fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2016 -DDAMASKVERSION=$DAMASKVERSION \
-qopenmp -qopenmp-threadprivate=compat\
-I$MARC_SOURCE/${BCS_DIR}/common -I$MARC_SOURCE/mumpssolver/include $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"

2
installation/mods_MarcMentat/installation.txt
View File

@ -16,7 +16,7 @@ The Intel Fortran compiler needs to be installed.
APPENDIX:
The structure of this directory should be (VERSION = 2010.2 or 2011 or 2012 or 2013 or 2014):
The structure of this directory should be (VERSION = 20XX or 20XX.Y)
./installation.txt this text
./apply_MPIE_modifications script file to apply modifications to the installation

0
lib/damask/util.py Executable file → Normal file
View File

1
src/CMakeLists.txt
View File

@ -74,6 +74,7 @@ add_library (PLASTIC OBJECT
"plastic_disloUCLA.f90"
"plastic_isotropic.f90"
"plastic_phenopowerlaw.f90"
"plastic_kinematichardening.f90"
"plastic_nonlocal.f90"
"plastic_none.f90")
add_dependencies(PLASTIC DAMASK_HELPERS)

3
src/CPFEM.f90
View File

@ -162,6 +162,7 @@ subroutine CPFEM_init
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
flush(6)
endif mainProcess
! initialize stress and jacobian to zero
@ -242,8 +243,8 @@ subroutine CPFEM_init
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
write(6,'(a32,1x,6(i8,1x),/)') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
write(6,'(a32,l1)') 'symmetricSolver: ', symmetricSolver
endif
flush(6)
endif
end subroutine CPFEM_init

51
src/DAMASK_spectral.f90
View File

@ -12,6 +12,8 @@ program DAMASK_spectral
compiler_version, &
compiler_options
#endif
#include <petsc/finclude/petscsys.h>
use PETScsys
use prec, only: &
pInt, &
pLongInt, &
@ -70,7 +72,6 @@ program DAMASK_spectral
DAMAGE_nonlocal_ID
use spectral_utilities, only: &
utilities_init, &
utilities_destroy, &
tSolutionState, &
tLoadCase, &
cutBack, &
@ -84,11 +85,8 @@ program DAMASK_spectral
use spectral_damage
use spectral_thermal
implicit none
#include <petsc/finclude/petscsys.h>
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
@ -143,18 +141,11 @@ program DAMASK_spectral
integer(pInt), parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742
integer(pInt), parameter :: maxRealOut = maxByteOut/pReal
integer(pLongInt), dimension(2) :: outputIndex
PetscErrorCode :: ierr
integer :: ierr
external :: &
quit, &
MPI_file_open, &
MPI_file_close, &
MPI_file_seek, &
MPI_file_get_position, &
MPI_file_write, &
MPI_abort, &
MPI_finalize, &
MPI_allreduce, &
PETScFinalize
quit
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
@ -442,10 +433,9 @@ program DAMASK_spectral
do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & ! QUESTION: why not starting i at 0 instead of murky 1?
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(resUnit, &
reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
[(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt), &
int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)), &
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
enddo
@ -635,7 +625,7 @@ program DAMASK_spectral
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),&
[(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt),&
int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)),&
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write')
enddo
@ -682,22 +672,12 @@ end program DAMASK_spectral
!> stderr. Exit code 3 signals no severe problems, but some increments did not converge
!--------------------------------------------------------------------------------------------------
subroutine quit(stop_id)
#include <petsc/finclude/petscsys.h>
use MPI
use prec, only: &
pInt
use spectral_mech_Basic, only: &
Basic_destroy
use spectral_mech_Polarisation, only: &
Polarisation_destroy
use spectral_damage, only: &
spectral_damage_destroy
use spectral_thermal, only: &
spectral_thermal_destroy
use spectral_utilities, only: &
utilities_destroy
implicit none
#include <petsc/finclude/petscsys.h>
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
integer(pInt) :: error = 0_pInt
@ -705,14 +685,7 @@ subroutine quit(stop_id)
logical :: ErrorInQuit
external :: &
PETScFinalize, &
MPI_finalize
call Basic_destroy()
call Polarisation_destroy()
call spectral_damage_destroy()
call spectral_thermal_destroy()
call utilities_destroy()
PETScFinalize
call PETScFinalize(ierr)
if (ierr /= 0) write(6,'(a)') ' Error in PETScFinalize'

1
src/commercialFEM_fileList.f90
View File

@ -28,6 +28,7 @@
#include "plastic_none.f90"
#include "plastic_isotropic.f90"
#include "plastic_phenopowerlaw.f90"
#include "plastic_kinematichardening.f90"
#include "plastic_dislotwin.f90"
#include "plastic_disloUCLA.f90"
#include "plastic_nonlocal.f90"

37
src/constitutive.f90
View File

@ -74,6 +74,7 @@ subroutine constitutive_init()
PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, &
PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_kinehardening_ID, &
PLASTICITY_dislotwin_ID, &
PLASTICITY_disloucla_ID, &
PLASTICITY_nonlocal_ID ,&
@ -95,6 +96,7 @@ subroutine constitutive_init()
PLASTICITY_NONE_label, &
PLASTICITY_ISOTROPIC_label, &
PLASTICITY_PHENOPOWERLAW_label, &
PLASTICITY_KINEHARDENING_label, &
PLASTICITY_DISLOTWIN_label, &
PLASTICITY_DISLOUCLA_label, &
PLASTICITY_NONLOCAL_label, &
@ -113,6 +115,7 @@ subroutine constitutive_init()
use plastic_none
use plastic_isotropic
use plastic_phenopowerlaw
use plastic_kinehardening
use plastic_dislotwin
use plastic_disloucla
use plastic_nonlocal
@ -156,6 +159,7 @@ subroutine constitutive_init()
if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then
@ -214,6 +218,11 @@ subroutine constitutive_init()
thisNoutput => plastic_phenopowerlaw_Noutput
thisOutput => plastic_phenopowerlaw_output
thisSize => plastic_phenopowerlaw_sizePostResult
case (PLASTICITY_KINEHARDENING_ID) plasticityType
outputName = PLASTICITY_KINEHARDENING_label
thisNoutput => plastic_kinehardening_Noutput
thisOutput => plastic_kinehardening_output
thisSize => plastic_kinehardening_sizePostResult
case (PLASTICITY_DISLOTWIN_ID) plasticityType
outputName = PLASTICITY_DISLOTWIN_label
thisNoutput => plastic_dislotwin_Noutput
@ -472,6 +481,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_KINEHARDENING_ID, &
PLASTICITY_DISLOTWIN_ID, &
PLASTICITY_DISLOUCLA_ID, &
PLASTICITY_NONLOCAL_ID
@ -479,6 +489,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
plastic_isotropic_LpAndItsTangent
use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_LpAndItsTangent
use plastic_kinehardening, only: &
plastic_kinehardening_LpAndItsTangent
use plastic_dislotwin, only: &
plastic_dislotwin_LpAndItsTangent
use plastic_disloucla, only: &
@ -525,6 +537,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, &
temperature(ho)%p(tme),ip,el)
@ -717,7 +731,7 @@ end function constitutive_initialFi
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic deformation gradient depending on the selected elastic law (so far no case switch
!! because only hooke is implemented
!! because only Hooke is implemented
!--------------------------------------------------------------------------------------------------
subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
use prec, only: &
@ -844,6 +858,7 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, &
PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_kinehardening_ID, &
PLASTICITY_dislotwin_ID, &
PLASTICITY_disloucla_ID, &
PLASTICITY_nonlocal_ID, &
@ -855,6 +870,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
plastic_isotropic_dotState
use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_dotState
use plastic_kinehardening, only: &
plastic_kinehardening_dotState
use plastic_dislotwin, only: &
plastic_dislotwin_dotState
use plastic_disloucla, only: &
@ -905,6 +922,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
call plastic_isotropic_dotState (Tstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_dotState(Tstar_v,ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dotState (Tstar_v,temperature(ho)%p(tme), &
ipc,ip,el)
@ -959,10 +978,13 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
phase_source, &
phase_Nsources, &
material_phase, &
PLASTICITY_KINEHARDENING_ID, &
PLASTICITY_NONLOCAL_ID, &
SOURCE_damage_isoBrittle_ID, &
SOURCE_vacancy_irradiation_ID, &
SOURCE_vacancy_thermalfluc_ID
use plastic_kinehardening, only: &
plastic_kinehardening_deltaState
use plastic_nonlocal, only: &
plastic_nonlocal_deltaState
use source_damage_isoBrittle, only: &
@ -991,9 +1013,12 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) &
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
if(phase_plasticity(material_phase(ipc,ip,el)) == PLASTICITY_NONLOCAL_ID) &
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(Tstar_v,ip,el)
end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
@ -1043,6 +1068,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_KINEHARDENING_ID, &
PLASTICITY_DISLOTWIN_ID, &
PLASTICITY_DISLOUCLA_ID, &
PLASTICITY_NONLOCAL_ID, &
@ -1054,6 +1080,8 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
plastic_isotropic_postResults
use plastic_phenopowerlaw, only: &
plastic_phenopowerlaw_postResults
use plastic_kinehardening, only: &
plastic_kinehardening_postResults
use plastic_dislotwin, only: &
plastic_dislotwin_postResults
use plastic_disloucla, only: &
@ -1102,6 +1130,9 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_kinehardening_postResults(Tstar_v,ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_dislotwin_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el)

82
src/crystallite.f90
View File

@ -986,7 +986,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst ! still on track or already done (beyond repair)
!$OMP FLUSH(crystallite_todo)
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt) then
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. c == debug_g) &
.or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) then
if (crystallite_todo(c,i,e)) then
write(6,'(a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent &
&with new crystallite_subStep: ',&
@ -1042,16 +1044,25 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif timeSyncing2
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then
write(6,'(/,a,e12.5)') '<< CRYST >> min(subStep) ',minval(crystallite_subStep)
write(6,'(a,e12.5)') '<< CRYST >> max(subStep) ',maxval(crystallite_subStep)
write(6,'(a,e12.5)') '<< CRYST >> min(subFrac) ',minval(crystallite_subFrac)
write(6,'(a,e12.5,/)') '<< CRYST >> max(subFrac) ',maxval(crystallite_subFrac)
write(6,'(/,a,f8.5)') '<< CRYST >> min(subStep) ',minval(crystallite_subStep)
write(6,'(a,f8.5)') '<< CRYST >> max(subStep) ',maxval(crystallite_subStep)
write(6,'(a,f8.5)') '<< CRYST >> min(subFrac) ',minval(crystallite_subFrac)
write(6,'(a,f8.5,/)') '<< CRYST >> max(subFrac) ',maxval(crystallite_subFrac)
flush(6)
if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt) then
write(6,'(/,a,f8.5,1x,a,1x,f8.5,1x,a)') '<< CRYST >> subFrac + subStep = ',&
crystallite_subFrac(debug_g,debug_i,debug_e),'+',crystallite_subStep(debug_g,debug_i,debug_e),'@selective'
flush(6)
endif
endif
! --- integrate --- requires fully defined state array (basic + dependent state)
if (any(crystallite_todo)) then
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then
write(6,'(/,a,i3)') '<< CRYST >> using state integrator ',numerics_integrator(numerics_integrationMode)
flush(6)
endif
select case(numerics_integrator(numerics_integrationMode))
case(1_pInt)
call crystallite_integrateStateFPI()
@ -2702,6 +2713,9 @@ subroutine crystallite_integrateStateFPI()
singleRun = (eIter(1) == eIter(2) .and. iIter(1,eIter(1)) == iIter(2,eIter(2)))
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo at start of state integration'
!--------------------------------------------------------------------------------------------------
! initialize dotState
if (.not. singleRun) then
@ -2754,6 +2768,8 @@ subroutine crystallite_integrateStateFPI()
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,*) '<< CRYST >> dotstate ',plasticState(p)%dotState(:,c)
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
@ -2767,6 +2783,9 @@ subroutine crystallite_integrateStateFPI()
!$OMP ENDDO
! --- UPDATE STATE ---
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after preguess of state'
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
@ -2822,6 +2841,9 @@ subroutine crystallite_integrateStateFPI()
! --- STRESS INTEGRATION ---
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo before stress integration'
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
@ -2976,7 +2998,11 @@ subroutine crystallite_integrateStateFPI()
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3,/)') '<< CRYST >> update state at el ip g ',e,i,g
write(6,'(a,f6.1,/)') '<< CRYST >> plasticstatedamper ',plasticStatedamper
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> plastic state residuum',plasticStateResiduum(1:mySizePlasticDotState)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> plastic state residuum',&
abs(plasticStateResiduum(1:mySizePlasticDotState))
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> abstol dotstate',plasticState(p)%aTolState(1:mySizePlasticDotState)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> reltol dotstate',rTol_crystalliteState* &
abs(tempPlasticState(1:mySizePlasticDotState))
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state',tempPlasticState(1:mySizePlasticDotState)
endif
#endif
@ -3036,7 +3062,7 @@ subroutine crystallite_integrateStateFPI()
!$OMP END PARALLEL
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
' grains converged after state integration #', NiterationState
@ -3152,8 +3178,8 @@ logical function crystallite_stateJump(ipc,ip,el)
write(6,'(a,i8,1x,i2,1x,i3, /)') '<< CRYST >> update state at el ip ipc ',el,ip,ipc
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', &
plasticState(p)%state(myOffsetSourceDeltaState + 1_pInt : &
myOffsetSourceDeltaState + mySizeSourceDeltaState,c)
plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : &
myOffsetPlasticDeltaState + mySizePlasticDeltaState,c)
endif
#endif
@ -3330,7 +3356,6 @@ logical function crystallite_integrateStress(&
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress at el ip ipc ',el,ip,ipc
#endif
!* only integrate over fraction of timestep?
if (present(timeFraction)) then
@ -3417,7 +3442,7 @@ logical function crystallite_integrateStress(&
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached loop limit',nStress, &
' at el (elFE) ip ipc ', el,mesh_element(1,el),ip,ipc
' at el (elFE) ip ipc ', el,'(',mesh_element(1,el),')',ip,ipc
#endif
return
endif loopsExeced
@ -3426,7 +3451,8 @@ logical function crystallite_integrateStress(&
B = math_I3 - dt*Lpguess
Fe = math_mul33x33(math_mul33x33(A,B), invFi_new) ! current elastic deformation tensor
call constitutive_TandItsTangent(Tstar, dT_dFe3333, dT_dFi3333, Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
call constitutive_TandItsTangent(Tstar, dT_dFe3333, dT_dFi3333, &
Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
Tstar_v = math_Mandel33to6(Tstar)
!* calculate plastic velocity gradient and its tangent from constitutive law
@ -3434,6 +3460,17 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i3,/)') '<< CRYST >> stress iteration ', NiterationStressLp
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Fi', math_transpose33(Fi_new)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Fe', math_transpose33(Fe)
write(6,'(a,/,6(e20.10,1x))') '<< CRYST >> Tstar', Tstar_v
endif
#endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT3333, dLp_dFi3333, &
Tstar_v, Fi_new, ipc, ip, el)
@ -3451,9 +3488,7 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i3,/)') '<< CRYST >> stress iteration ', NiterationStressLp
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
endif
#endif
@ -3483,6 +3518,13 @@ logical function crystallite_integrateStress(&
else ! not converged and residuum not improved...
steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplengthLp * deltaLp
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,1x,f7.4)') '<< CRYST >> linear search for Lpguess with step', steplengthLp
endif
#endif
cycle LpLoop
endif
@ -3496,6 +3538,16 @@ logical function crystallite_integrateStress(&
dFe_dLp3333 = - dt * dFe_dLp3333
dRLp_dLp = math_identity2nd(9_pInt) &
- math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333))
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dLp_dT', math_Plain3333to99(dLp_dT3333)
write(6,'(a,1x,e20.10)') '<< CRYST >> dLp_dT norm', norm2(math_Plain3333to99(dLp_dT3333))
write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dRLp_dLp', dRLp_dLp - math_identity2nd(9_pInt)
write(6,'(a,1x,e20.10)') '<< CRYST >> dRLp_dLp norm', norm2(dRLp_dLp - math_identity2nd(9_pInt))
endif
#endif
dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine
work = math_plain33to9(residuumLp)
call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp

14
src/debug.f90
View File

@ -443,13 +443,15 @@ subroutine debug_info
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
endif debugOutputHomog
debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0) then
write(6,'(2/,a,/)') ' Extreme values of returned stress and jacobian'
debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 &
.and. any(debug_stressMinLocation /= 0_pInt) &
.and. any(debug_stressMaxLocation /= 0_pInt) ) then
write(6,'(2/,a,/)') ' Extreme values of returned stress and Jacobian'
write(6,'(a39)') ' value el ip'
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' Jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
endif debugOutputCPFEM
!$OMP END CRITICAL (write2out)

0
src/homogenization_RGC.f90 Executable file → Normal file
View File

5
src/material.f90
View File

@ -25,6 +25,7 @@ module material
PLASTICITY_none_label = 'none', &
PLASTICITY_isotropic_label = 'isotropic', &
PLASTICITY_phenopowerlaw_label = 'phenopowerlaw', &
PLASTICITY_kinehardening_label = 'kinehardening', &
PLASTICITY_dislotwin_label = 'dislotwin', &
PLASTICITY_disloucla_label = 'disloucla', &
PLASTICITY_nonlocal_label = 'nonlocal', &
@ -72,6 +73,7 @@ module material
PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, &
PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_kinehardening_ID, &
PLASTICITY_dislotwin_ID, &
PLASTICITY_disloucla_ID, &
PLASTICITY_nonlocal_ID
@ -308,6 +310,7 @@ module material
PLASTICITY_none_ID, &
PLASTICITY_isotropic_ID, &
PLASTICITY_phenopowerlaw_ID, &
PLASTICITY_kinehardening_ID, &
PLASTICITY_dislotwin_ID, &
PLASTICITY_disloucla_ID, &
PLASTICITY_nonlocal_ID, &
@ -983,6 +986,8 @@ subroutine material_parsePhase(fileUnit,myPart)
phase_plasticity(section) = PLASTICITY_ISOTROPIC_ID
case (PLASTICITY_PHENOPOWERLAW_label)
phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID
case (PLASTICITY_KINEHARDENING_label)
phase_plasticity(section) = PLASTICITY_KINEHARDENING_ID
case (PLASTICITY_DISLOTWIN_label)
phase_plasticity(section) = PLASTICITY_DISLOTWIN_ID
case (PLASTICITY_DISLOUCLA_label)

557
src/math.f90
View File

@ -74,6 +74,7 @@ module math
public :: &
math_init, &
math_qsort, &
math_expand, &
math_range, &
math_identity2nd, &
math_identity4th, &
@ -161,10 +162,8 @@ module math
math_rotate_forward3333, &
math_limit
private :: &
halton, &
halton_memory, &
halton_ndim_set, &
halton_seed_set
math_check, &
halton
contains
@ -216,10 +215,6 @@ subroutine math_init
write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest
call random_seed(put = randInit)
call halton_seed_set(int(randInit(1), pInt))
call halton_ndim_set(3_pInt)
call math_check()
end subroutine math_init
@ -283,7 +278,7 @@ subroutine math_check
endif
! +++ check rotation sense of q and R +++
call halton(3_pInt,v) ! random vector
v = halton([2_pInt,8_pInt,5_pInt]) ! random vector
R = math_qToR(q)
if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then
write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v'
@ -1237,7 +1232,7 @@ function math_qRand()
real(pReal), dimension(4) :: math_qRand
real(pReal), dimension(3) :: rnd
call halton(3_pInt,rnd)
rnd = halton([8_pInt,4_pInt,9_pInt])
math_qRand = [cos(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)), &
sin(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
cos(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
@ -1760,7 +1755,7 @@ function math_sampleRandomOri()
implicit none
real(pReal), dimension(3) :: math_sampleRandomOri, rnd
call halton(3_pInt,rnd)
rnd = halton([1_pInt,7_pInt,3_pInt])
math_sampleRandomOri = [rnd(1)*2.0_pReal*PI, &
acos(2.0_pReal*rnd(2)-1.0_pReal), &
rnd(3)*2.0_pReal*PI]
@ -1769,118 +1764,96 @@ end function math_sampleRandomOri
!--------------------------------------------------------------------------------------------------
!> @brief draw a random sample from Gauss component with noise (in radians) half-width
!> @brief draw a sample from an Gaussian distribution around given orientation and Full Width
! at Half Maximum (FWHM)
!> @details: A uniform misorientation (limited to 2*FWHM) is sampled followed by convolution with
! a Gausian distribution
!--------------------------------------------------------------------------------------------------
function math_sampleGaussOri(center,noise)
use prec, only: &
tol_math_check
function math_sampleGaussOri(center,FWHM)
implicit none
real(pReal), intent(in) :: noise
real(pReal), intent(in) :: FWHM
real(pReal), dimension(3), intent(in) :: center
real(pReal) :: cosScatter,scatter
real(pReal), dimension(3) :: math_sampleGaussOri, disturb
real(pReal), dimension(3), parameter :: ORIGIN = 0.0_pReal
real(pReal), dimension(5) :: rnd
real(pReal) :: angle
real(pReal), dimension(3) :: math_sampleGaussOri, axis
real(pReal), dimension(4) :: rnd
real(pReal), dimension(3,3) :: R
noScatter: if (abs(noise) < tol_math_check) then
if (FWHM < 0.1_pReal*INRAD) then
math_sampleGaussOri = center
else noScatter
! Helming uses different distribution with Bessel functions
! therefore the gauss scatter width has to be scaled differently
scatter = 0.95_pReal * noise
cosScatter = cos(scatter)
else
GaussConvolution: do
rnd = halton([8_pInt,3_pInt,6_pInt,11_pInt])
axis(1) = rnd(1)*2.0_pReal-1.0_pReal ! uniform on [-1,1]
axis(2:3) = [sqrt(1.0-axis(1)**2.0_pReal)*cos(rnd(2)*2.0*PI),&
sqrt(1.0-axis(1)**2.0_pReal)*sin(rnd(2)*2.0*PI)] ! random axis
angle = (rnd(3)-0.5_pReal)*4.0_pReal*FWHM ! rotation by [0, +-2 FWHM]
R = math_axisAngleToR(axis,angle)
angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R))
if (rnd(4) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) exit ! rejection sampling (Gaussian)
enddo GaussConvolution
math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center)))
endif
do
call halton(5_pInt,rnd)
rnd(1:3) = 2.0_pReal*rnd(1:3)-1.0_pReal ! expand 1:3 to range [-1,+1]
disturb = [ scatter * rnd(1), & ! phi1
sign(1.0_pReal,rnd(2))*acos(cosScatter+(1.0_pReal-cosScatter)*rnd(4)), & ! Phi
scatter * rnd(3)] ! phi2
if (rnd(5) <= exp(-1.0_pReal*(math_EulerMisorientation(ORIGIN,disturb)/scatter)**2_pReal)) exit
enddo
math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center)))
endif noScatter
end function math_sampleGaussOri
!--------------------------------------------------------------------------------------------------
!> @brief draw a random sample from Fiber component with noise (in radians)
!--------------------------------------------------------------------------------------------------
function math_sampleFiberOri(alpha,beta,noise)
use prec, only: &
tol_math_check
!> @brief draw a sample from an Gaussian distribution around given fiber texture and Full Width
! at Half Maximum (FWHM)
!-------------------------------------------------------------------------------------------------
function math_sampleFiberOri(alpha,beta,FWHM)
implicit none
real(pReal), dimension(3) :: math_sampleFiberOri, fiberInC,fiberInS,axis
real(pReal), dimension(2), intent(in) :: alpha,beta
real(pReal), dimension(6) :: rnd
real(pReal), dimension(3,3) :: oRot,fRot,pRot
real(pReal) :: noise, scatter, cos2Scatter, angle
integer(pInt), dimension(2,3), parameter :: ROTMAP = reshape([2_pInt,3_pInt,&
3_pInt,1_pInt,&
1_pInt,2_pInt],[2,3])
integer(pInt) :: i
real(pReal), intent(in) :: FWHM
real(pReal), dimension(3) :: math_sampleFiberOri, &
fInC,& !< fiber axis in crystal coordinate system
fInS,& !< fiber axis in sample coordinate system
u
real(pReal), dimension(3) :: rnd
real(pReal), dimension(:),allocatable :: a !< 2D vector to tilt
integer(pInt), dimension(:),allocatable :: idx !< components of 2D vector
real(pReal), dimension(3,3) :: R !< Rotation matrix (composed of three components)
real(pReal):: angle,c
integer(pInt):: j,& !< index of smallest component
i
! Helming uses different distribution with Bessel functions
! therefore the gauss scatter width has to be scaled differently
scatter = 0.95_pReal * noise
cos2Scatter = cos(2.0_pReal*scatter)
fInC = [sin(alpha(1))*cos(alpha(2)), sin(alpha(1))*sin(alpha(2)), cos(alpha(1))]
fInS = [sin(beta(1))*cos(beta(2)), sin(beta(1))*sin(beta(2)), cos(beta(1))]
! fiber axis in crystal coordinate system
fiberInC = [ sin(alpha(1))*cos(alpha(2)) , &
sin(alpha(1))*sin(alpha(2)), &
cos(alpha(1))]
! fiber axis in sample coordinate system
fiberInS = [ sin(beta(1))*cos(beta(2)), &
sin(beta(1))*sin(beta(2)), &
cos(beta(1))]
R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system
! ---# rotation matrix from sample to crystal system #---
angle = -acos(dot_product(fiberInC,fiberInS))
if(abs(angle) > tol_math_check) then
! rotation axis between sample and crystal system (cross product)
forall(i=1_pInt:3_pInt) axis(i) = fiberInC(ROTMAP(1,i))*fiberInS(ROTMAP(2,i))-fiberInC(ROTMAP(2,i))*fiberInS(ROTMAP(1,i))
oRot = math_EulerAxisAngleToR(math_crossproduct(fiberInC,fiberInS),angle)
rnd = halton([7_pInt,10_pInt,3_pInt])
R = math_mul33x33(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis
if (FWHM > 0.1_pReal*INRAD) then
reducedTo2D: do i=1_pInt,3_pInt
if (i /= minloc(abs(fInS),1)) then
a=[a,fInS(i)]
idx=[idx,i]
else
oRot = math_I3
j = i
endif
enddo reducedTo2D
GaussConvolution: do
angle = (rnd(2)-0.5_pReal)*4.0_pReal*FWHM ! rotation by [0, +-2 FWHM]
! solve cos(angle) = dot_product(fInS,u) under the assumption that their smallest component is the same
c = cos(angle)-fInS(j)**2
u(idx(2)) = -(2.0_pReal*c*a(2) + sqrt(4*((c*a(2))**2-sum(a**2)*(c**2-a(1)**2*(1-fInS(j)**2)))))/&
(2*sum(a**2))
u(idx(1)) = sqrt(1-u(idx(2))**2-fInS(j)**2)
u(j) = fInS(j)
! ---# rotation matrix about fiber axis (random angle) #---
do
call halton(6_pInt,rnd)
fRot = math_EulerAxisAngleToR(fiberInS,rnd(1)*2.0_pReal*pi)
! ---# rotation about random axis perpend to fiber #---
! random axis pependicular to fiber axis
axis(1:2) = rnd(2:3)
if (abs(fiberInS(3)) > tol_math_check) then
axis(3)=-(axis(1)*fiberInS(1)+axis(2)*fiberInS(2))/fiberInS(3)
else if(abs(fiberInS(2)) > tol_math_check) then
axis(3)=axis(2)
axis(2)=-(axis(1)*fiberInS(1)+axis(3)*fiberInS(3))/fiberInS(2)
else if(abs(fiberInS(1)) > tol_math_check) then
axis(3)=axis(1)
axis(1)=-(axis(2)*fiberInS(2)+axis(3)*fiberInS(3))/fiberInS(1)
end if
! scattered rotation angle
if (noise > 0.0_pReal) then
angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(4))
if (rnd(5) <= exp(-1.0_pReal*(angle/scatter)**2.0_pReal)) exit
else
angle = 0.0_pReal
rejectionSampling: if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) then
R = math_mul33x33(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component
exit
endif rejectionSampling
rnd = halton([7_pInt,10_pInt,3_pInt])
enddo GaussConvolution
endif
enddo
if (rnd(6) <= 0.5) angle = -angle
pRot = math_EulerAxisAngleToR(axis,angle)
! ---# apply the three rotations #---
math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot)))
math_sampleFiberOri = math_RtoEuler(R)
end function math_sampleFiberOri
@ -1902,18 +1875,17 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width)
if (abs(stddev) < tol_math_check) then
math_sampleGaussVar = meanvalue
return
endif
else
myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given
do
call halton(2_pInt, rnd)
rnd = halton([6_pInt,2_pInt])
scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal)
if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn
enddo
math_sampleGaussVar = scatter * stddev
endif
end function math_sampleGaussVar
@ -2284,172 +2256,36 @@ pure function math_invariantsSym33(m)
end function math_invariantsSym33
!--------------------------------------------------------------------------------------------------
!> @brief computes the next element in the Halton sequence.
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
subroutine halton(ndim, r)
implicit none
integer(pInt), intent(in) :: ndim !< dimension of the element
real(pReal), intent(out), dimension(ndim) :: r !< next element of the current Halton sequence
integer(pInt), dimension(ndim) :: base
integer(pInt) :: seed
integer(pInt), dimension(1) :: value_halton
call halton_memory ('GET', 'SEED', 1_pInt, value_halton)
seed = value_halton(1)
call halton_memory ('GET', 'BASE', ndim, base)
call i_to_halton (seed, base, ndim, r)
value_halton(1) = 1_pInt
call halton_memory ('INC', 'SEED', 1_pInt, value_halton)
!--------------------------------------------------------------------------------------------------
contains
!-------------------------------------------------------------------------------------------------
!> @brief computes an element of a Halton sequence.
!> @details Only the absolute value of SEED is considered. SEED = 0 is allowed, and returns R = 0.
!> @details Halton Bases should be distinct prime numbers. This routine only checks that each base
!> @details is greater than 1.
!> @author John Burkardt
!> @author Martin Diehl
!> @details Incrementally increasing elements of the Halton sequence for given bases (> 0)
!> @details Reference:
!> @details J.H. Halton: On the efficiency of certain quasi-random sequences of points in evaluating
!> @details multi-dimensional integrals, Numerische Mathematik, Volume 2, pages 84-90, 1960.
!> @author John Burkardt
!-------------------------------------------------------------------------------------------------
subroutine i_to_halton (seed, base, ndim, r)
use IO, only: &
IO_error
implicit none
integer(pInt), intent(in) :: &
ndim, & !< dimension of the sequence
seed !< index of the desired element
integer(pInt), intent(in), dimension(ndim) :: base !< Halton bases
real(pReal), intent(out), dimension(ndim) :: r !< the SEED-th element of the Halton sequence for the given bases
real(pReal), dimension(ndim) :: base_inv
integer(pInt), dimension(ndim) :: &
digit, &
seed2
seed2 = abs(seed)
r = 0.0_pReal
if (any (base(1:ndim) <= 1_pInt)) call IO_error(error_ID=405_pInt)
base_inv(1:ndim) = 1.0_pReal / real (base(1:ndim), pReal)
do while ( any ( seed2(1:ndim) /= 0_pInt) )
digit(1:ndim) = mod ( seed2(1:ndim), base(1:ndim))
r(1:ndim) = r(1:ndim) + real ( digit(1:ndim), pReal) * base_inv(1:ndim)
base_inv(1:ndim) = base_inv(1:ndim) / real ( base(1:ndim), pReal)
seed2(1:ndim) = seed2(1:ndim) / base(1:ndim)
enddo
end subroutine i_to_halton
end subroutine halton
!--------------------------------------------------------------------------------------------------
!> @brief sets or returns quantities associated with the Halton sequence.
!> @details If action_halton is 'SET' and action_halton is 'BASE', then NDIM is input, and
!> @details is the number of entries in value_halton to be put into BASE.
!> @details If action_halton is 'SET', then on input, value_halton contains values to be assigned
!> @details to the internal variable.
!> @details If action_halton is 'GET', then on output, value_halton contains the values of
!> @details the specified internal variable.
!> @details If action_halton is 'INC', then on input, value_halton contains the increment to
!> @details be added to the specified internal variable.
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
use IO, only: &
IO_lc
implicit none
character(len = *), intent(in) :: &
action_halton, & !< desired action: GET the value of a particular quantity, SET the value of a particular quantity, INC the value of a particular quantity (only for SEED)
name_halton !< name of the quantity: BASE: Halton base(s), NDIM: spatial dimension, SEED: current Halton seed
integer(pInt), dimension(*), intent(inout) :: value_halton
integer(pInt), allocatable, save, dimension(:) :: base
logical, save :: first_call = .true.
integer(pInt), intent(in) :: ndim !< dimension of the quantity
integer(pInt), save :: ndim_save = 0_pInt, seed = 1_pInt
integer(pInt) :: i
if (first_call) then
ndim_save = 1_pInt
allocate(base(ndim_save))
base(1) = 2_pInt
first_call = .false.
endif
!--------------------------------------------------------------------------------------------------
! Set
actionHalton: if(IO_lc(action_halton(1:1)) == 's') then
nameSet: if(IO_lc(name_halton(1:1)) == 'b') then
if(ndim_save /= ndim) ndim_save = ndim
base = value_halton(1:ndim)
elseif(IO_lc(name_halton(1:1)) == 'n') then nameSet
if(ndim_save /= value_halton(1)) then
ndim_save = value_halton(1)
base = [(prime(i),i=1_pInt,ndim_save)]
else
ndim_save = value_halton(1)
endif
elseif(IO_lc(name_halton(1:1)) == 's') then nameSet
seed = value_halton(1)
endif nameSet
!--------------------------------------------------------------------------------------------------
! Get
elseif(IO_lc(action_halton(1:1)) == 'g') then actionHalton
nameGet: if(IO_lc(name_halton(1:1)) == 'b') then
if(ndim /= ndim_save) then
ndim_save = ndim
base = [(prime(i),i=1_pInt,ndim_save)]
endif
value_halton(1:ndim_save) = base(1:ndim_save)
elseif(IO_lc(name_halton(1:1)) == 'n') then nameGet
value_halton(1) = ndim_save
elseif(IO_lc(name_halton(1:1)) == 's') then nameGet
value_halton(1) = seed
endif nameGet
!--------------------------------------------------------------------------------------------------
! Increment
elseif(IO_lc(action_halton(1:1)) == 'i') then actionHalton
if(IO_lc(name_halton(1:1)) == 's') seed = seed + value_halton(1)
endif actionHalton
!--------------------------------------------------------------------------------------------------
contains
!--------------------------------------------------------------------------------------------------
!> @brief returns any of the first 1500 prime numbers.
!> @details n = 0 is legal, returning PRIME = 1.
!> @details Reference:
!> @details Reference for prime numbers:
!> @details Milton Abramowitz and Irene Stegun: Handbook of Mathematical Functions,
!> @details US Department of Commerce, 1964, pages 870-873.
!> @details Daniel Zwillinger: CRC Standard Mathematical Tables and Formulae,
!> @details 30th Edition, CRC Press, 1996, pages 95-98.
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
integer(pInt) function prime(n)
use IO, only: &
IO_error
!-------------------------------------------------------------------------------------------------
function halton(bases)
implicit none
integer(pInt), intent(in) :: n !< index of the desired prime number
integer(pInt), dimension(0:1500), parameter :: &
npvec = int([&
integer(pInt), intent(in), dimension(:):: &
bases !< bases (prime number ID)
real(pReal), dimension(size(bases)) :: &
halton
integer(pInt), save :: &
current = 1_pInt
real(pReal), dimension(size(bases)) :: &
base_inv
integer(pInt), dimension(size(bases)) :: &
base, &
t
integer(pInt), dimension(0:1600), parameter :: &
prime = int([&
1, &
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, &
31, 37, 41, 43, 47, 53, 59, 61, 67, 71, &
@ -2596,7 +2432,7 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
! 1301:1400
10663, 10667, 10687, 10691, 10709, 10711, 10723, 10729, 10733, 10739, &
10753, 10771, 10781, 10789, 10799, 10831, 10837, 10847, 10853, 10859, &
10861, 10867, 10883, 10889, 10891, 10903, 10909, 19037, 10939, 10949, &
10861, 10867, 10883, 10889, 10891, 10903, 10909, 10937, 10939, 10949, &
10957, 10973, 10979, 10987, 10993, 11003, 11027, 11047, 11057, 11059, &
11069, 11071, 11083, 11087, 11093, 11113, 11117, 11119, 11131, 11149, &
11159, 11161, 11171, 11173, 11177, 11197, 11213, 11239, 11243, 11251, &
@ -2614,58 +2450,34 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
12227, 12239, 12241, 12251, 12253, 12263, 12269, 12277, 12281, 12289, &
12301, 12323, 12329, 12343, 12347, 12373, 12377, 12379, 12391, 12401, &
12409, 12413, 12421, 12433, 12437, 12451, 12457, 12473, 12479, 12487, &
12491, 12497, 12503, 12511, 12517, 12527, 12539, 12541, 12547, 12553],pInt)
12491, 12497, 12503, 12511, 12517, 12527, 12539, 12541, 12547, 12553, &
! 1501:1600
12569, 12577, 12583, 12589, 12601, 12611, 12613, 12619, 12637, 12641, &
12647, 12653, 12659, 12671, 12689, 12697, 12703, 12713, 12721, 12739, &
12743, 12757, 12763, 12781, 12791, 12799, 12809, 12821, 12823, 12829, &
12841, 12853, 12889, 12893, 12899, 12907, 12911, 12917, 12919, 12923, &
12941, 12953, 12959, 12967, 12973, 12979, 12983, 13001, 13003, 13007, &
13009, 13033, 13037, 13043, 13049, 13063, 13093, 13099, 13103, 13109, &
13121, 13127, 13147, 13151, 13159, 13163, 13171, 13177, 13183, 13187, &
13217, 13219, 13229, 13241, 13249, 13259, 13267, 13291, 13297, 13309, &
13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399, 13411, &
13417, 13421, 13441, 13451, 13457, 13463, 13469, 13477, 13487, 13499],pInt)
if (n < size(npvec)) then
prime = npvec(n)
else
call IO_error(error_ID=406_pInt)
end if
current = current + 1_pInt
end function prime
base = prime(bases)
base_inv = 1.0_pReal/real(base,pReal)
end subroutine halton_memory
halton = 0.0_pReal
t = current
do while (any( t /= 0_pInt) )
halton = halton + real(mod(t,base), pReal) * base_inv
base_inv = base_inv / real(base, pReal)
t = t / base
enddo
!--------------------------------------------------------------------------------------------------
!> @brief sets the dimension for a Halton sequence
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
subroutine halton_ndim_set(ndim)
implicit none
integer(pInt), intent(in) :: ndim !< dimension of the Halton vectors
integer(pInt) :: value_halton(1)
value_halton(1) = ndim
call halton_memory ('SET', 'NDIM', 1_pInt, value_halton)
end subroutine halton_ndim_set
!--------------------------------------------------------------------------------------------------
!> @brief sets the seed for the Halton sequence.
!> @details Calling HALTON repeatedly returns the elements of the Halton sequence in order,
!> @details starting with element number 1.
!> @details An internal counter, called SEED, keeps track of the next element to return. Each time
!> @details is computed, and then SEED is incremented by 1.
!> @details To restart the Halton sequence, it is only necessary to reset SEED to 1. It might also
!> @details be desirable to reset SEED to some other value. This routine allows the user to specify
!> @details any value of SEED.
!> @details The default value of SEED is 1, which restarts the Halton sequence.
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
subroutine halton_seed_set(seed)
implicit none
integer(pInt), parameter :: NDIM = 1_pInt
integer(pInt), intent(in) :: seed !< seed for the Halton sequence.
integer(pInt) :: value_halton(ndim)
value_halton(1) = seed
call halton_memory ('SET', 'SEED', NDIM, value_halton)
end subroutine halton_seed_set
end function halton
!--------------------------------------------------------------------------------------------------
@ -2820,4 +2632,135 @@ real(pReal) pure function math_limit(a, left, right)
end function math_limit
!--------------------------------------------------------------------------------------------------
!> @brief Modified Bessel I function of order 0
!> @author John Burkardt
!> @details original version available on https://people.sc.fsu.edu/~jburkardt/f_src/toms715/toms715.html
!--------------------------------------------------------------------------------------------------
real(pReal) function bessel_i0 (x)
use, intrinsic :: IEEE_ARITHMETIC
implicit none
real(pReal), intent(in) :: x
integer(pInt) :: i
real(pReal) :: sump_p, sump_q, xAbs, xx
real(pReal), parameter, dimension(15) :: p_small = real( &
[-5.2487866627945699800e-18, -1.5982226675653184646e-14, -2.6843448573468483278e-11, &
-3.0517226450451067446e-08, -2.5172644670688975051e-05, -1.5453977791786851041e-02, &
-7.0935347449210549190e+00, -2.4125195876041896775e+03, -5.9545626019847898221e+05, &
-1.0313066708737980747e+08, -1.1912746104985237192e+10, -8.4925101247114157499e+11, &
-3.2940087627407749166e+13, -5.5050369673018427753e+14, -2.2335582639474375249e+15], pReal)
real(pReal), parameter, dimension(5) :: q_small = real( &
[-3.7277560179962773046e+03, 6.5158506418655165707e+06, -6.5626560740833869295e+09, &
3.7604188704092954661e+12, -9.7087946179594019126e+14], pReal)
real(pReal), parameter, dimension(8) :: p_large = real( &
[-3.9843750000000000000e-01, 2.9205384596336793945e+00, -2.4708469169133954315e+00, &
4.7914889422856814203e-01, -3.7384991926068969150e-03, -2.6801520353328635310e-03, &
9.9168777670983678974e-05, -2.1877128189032726730e-06], pReal)
real(pReal), parameter, dimension(7) :: q_large = real( &
[-3.1446690275135491500e+01, 8.5539563258012929600e+01, -6.0228002066743340583e+01, &
1.3982595353892851542e+01, -1.1151759188741312645e+00, 3.2547697594819615062e-02, &
-5.5194330231005480228e-04], pReal)
xAbs = abs(x)
argRange: if (xAbs < 5.55e-17_pReal) then
bessel_i0 = 1.0_pReal
else if (xAbs < 15.0_pReal) then argRange
xx = xAbs**2.0_pReal
sump_p = p_small(1)
do i = 2, 15
sump_p = sump_p * xx + p_small(i)
end do
xx = xx - 225.0_pReal
sump_q = ((((xx+q_small(1))*xx+q_small(2))*xx+q_small(3))*xx+q_small(4))*xx+q_small(5)
bessel_i0 = sump_p / sump_q
else if (xAbs <= 713.986_pReal) then argRange
xx = 1.0_pReal / xAbs - 2.0_pReal/30.0_pReal
sump_p = ((((((p_large(1)*xx+p_large(2))*xx+p_large(3))*xx+p_large(4))*xx+ &
p_large(5))*xx+p_large(6))*xx+p_large(7))*xx+p_large(8)
sump_q = ((((((xx+q_large(1))*xx+q_large(2))*xx+q_large(3))*xx+ &
q_large(4))*xx+q_large(5))*xx+q_large(6))*xx+q_large(7)
bessel_i0 = sump_p / sump_q
avoidOverflow: if (xAbs > 698.986_pReal) then
bessel_i0 = ((bessel_i0*exp(xAbs-40.0_pReal)-p_large(1)*exp(xAbs-40.0_pReal))/sqrt(xAbs))*exp(40.0)
else avoidOverflow
bessel_i0 = ((bessel_i0*exp(xAbs)-p_large(1)*exp(xAbs))/sqrt(xAbs))
endif avoidOverflow
else argRange
bessel_i0 = IEEE_value(bessel_i0,IEEE_positive_inf)
end if argRange
end function bessel_i0
!--------------------------------------------------------------------------------------------------
!> @brief Modified Bessel I function of order 1
!> @author John Burkardt
!> @details original version available on https://people.sc.fsu.edu/~jburkardt/f_src/toms715/toms715.html
!--------------------------------------------------------------------------------------------------
real(pReal) function bessel_i1 (x)
use, intrinsic :: IEEE_ARITHMETIC
implicit none
real(pReal), intent(in) :: x
integer(pInt) :: i
real(pReal) :: sump_p, sump_q, xAbs, xx
real(pReal), dimension(15), parameter :: p_small = real( &
[-1.9705291802535139930e-19, -6.5245515583151902910e-16, -1.1928788903603238754e-12, &
-1.4831904935994647675e-09, -1.3466829827635152875e-06, -9.1746443287817501309e-04, &
-4.7207090827310162436e-01, -1.8225946631657315931e+02, -5.1894091982308017540e+04, &
-1.0588550724769347106e+07, -1.4828267606612366099e+09, -1.3357437682275493024e+11, &
-6.9876779648010090070e+12, -1.7732037840791591320e+14, -1.4577180278143463643e+15], pReal)
real(pReal), dimension(5), parameter :: q_small = real( &
[-4.0076864679904189921e+03, 7.4810580356655069138e+06, -8.0059518998619764991e+09, &
4.8544714258273622913e+12, -1.3218168307321442305e+15], pReal)
real(pReal), dimension(8), parameter :: p_large = real( &
[-6.0437159056137600000e-02, 4.5748122901933459000e-01, -4.2843766903304806403e-01, &
9.7356000150886612134e-02, -3.2457723974465568321e-03, -3.6395264712121795296e-04, &
1.6258661867440836395e-05, -3.6347578404608223492e-07], pReal)
real(pReal), dimension(6), parameter :: q_large = real( &
[-3.8806586721556593450e+00, 3.2593714889036996297e+00, -8.5017476463217924408e-01, &
7.4212010813186530069e-02, -2.2835624489492512649e-03, 3.7510433111922824643e-05], pReal)
real(pReal), parameter :: pbar = 3.98437500e-01
xAbs = abs(x)
argRange: if (xAbs < 5.55e-17_pReal) then
bessel_i1 = 0.5_pReal * xAbs
else if (xAbs < 15.0_pReal) then argRange
xx = xAbs**2.0_pReal
sump_p = p_small(1)
do i = 2, 15
sump_p = sump_p * xx + p_small(i)
end do
xx = xx - 225.0_pReal
sump_q = ((((xx+q_small(1))*xx+q_small(2))*xx+q_small(3))*xx+q_small(4)) * xx + q_small(5)
bessel_i1 = (sump_p / sump_q) * xAbs
else if (xAbs <= 713.986_pReal) then argRange
xx = 1.0_pReal / xAbs - 2.0_pReal/30.0_pReal
sump_p = ((((((p_large(1)*xx+p_large(2))*xx+p_large(3))*xx+p_large(4))*xx+&
p_large(5))*xx+p_large(6))*xx+p_large(7))*xx+p_large(8)
sump_q = (((((xx+q_large(1))*xx+q_large(2))*xx+q_large(3))*xx+ q_large(4))*xx+q_large(5))*xx+q_large(6)
bessel_i1 = sump_p / sump_q
avoidOverflow: if (xAbs > 698.986_pReal) then
bessel_i1 = ((bessel_i1 * exp(xAbs-40.0_pReal) + pbar * exp(xAbs-40.0_pReal)) / sqrt(xAbs)) * exp(40.0_pReal)
else avoidOverflow
bessel_i1 = ((bessel_i1 * exp(xAbs) + pbar * exp(xAbs)) / sqrt(xAbs))
endif avoidOverflow
else argRange
bessel_i1 = IEEE_value(bessel_i1,IEEE_positive_inf)
end if argRange
if (x < 0.0_pReal) bessel_i1 = -bessel_i1
end function bessel_i1
end module math

12
src/mesh.f90
View File

@ -118,11 +118,6 @@ module mesh
logical, private :: noPart !< for cases where the ABAQUS input file does not use part/assembly information
#endif
#ifdef Spectral
#include <petsc/finclude/petscsys.h>
include 'fftw3-mpi.f03'
#endif
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
! Hence, I suggest to prefix with "FE_"
@ -481,6 +476,10 @@ subroutine mesh_init(ip,el)
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
#ifdef Spectral
#include <petsc/finclude/petscsys.h>
use PETScsys
#endif
use DAMASK_interface
use IO, only: &
@ -516,6 +515,7 @@ subroutine mesh_init(ip,el)
implicit none
#ifdef Spectral
include 'fftw3-mpi.f03'
integer(C_INTPTR_T) :: devNull, local_K, local_K_offset
integer :: ierr, worldsize
#endif
@ -524,8 +524,6 @@ subroutine mesh_init(ip,el)
integer(pInt) :: j
logical :: myDebug
external :: MPI_comm_size
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"

11
src/numerics.f90
View File

@ -10,9 +10,6 @@ module numerics
implicit none
private
#ifdef PETSc
#include <petsc/finclude/petsc.h90>
#endif
character(len=64), parameter, private :: &
numerics_CONFIGFILE = 'numerics.config' !< name of configuration file
@ -216,6 +213,10 @@ subroutine numerics_init
IO_warning, &
IO_timeStamp, &
IO_EOF
#ifdef PETSc
#include <petsc/finclude/petscsys.h>
use petscsys
#endif
#if defined(Spectral) || defined(FEM)
!$ use OMP_LIB, only: omp_set_num_threads ! Use the standard conforming module file for omp if using the spectral solver
implicit none
@ -232,9 +233,7 @@ subroutine numerics_init
line
!$ character(len=6) DAMASK_NumThreadsString ! environment variable DAMASK_NUM_THREADS
external :: &
MPI_Comm_rank, &
MPI_Comm_size, &
MPI_Abort
PETScErrorF ! is called in the CHKERRQ macro
#ifdef PETSc
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)

2
src/plastic_isotropic.f90
View File

@ -62,7 +62,7 @@ module plastic_isotropic
end type
type, private :: tIsotropicAbsTol !< internal alias for abs tolerance in state
real(pReal), pointer :: & ! scalars along NipcMyInstance
real(pReal), pointer :: & ! scalars
flowstress, &
accumulatedShear
end type

969
src/plastic_kinematichardening.f90 Normal file
View File

@ -0,0 +1,969 @@
!--------------------------------------------------------------------------------------------------
!> @author Philip Eisenlohr, Michigan State University
!> @author Zhuowen Zhao, Michigan State University
!> @brief Introducing Voce-type kinematic hardening rule into crystal plasticity
!! formulation using a power law fitting
!--------------------------------------------------------------------------------------------------
module plastic_kinehardening
use prec, only: &
pReal,&
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_kinehardening_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: &
plastic_kinehardening_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_kinehardening_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
plastic_kinehardening_Noutput !< number of outputs per instance
integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_kinehardening_totalNslip !< no. of slip system used in simulation
integer(pInt), dimension(:,:), allocatable, private :: &
plastic_kinehardening_Nslip !< active number of slip systems per family (input parameter, per family)
enum, bind(c)
enumerator :: undefined_ID, &
crss_ID, & !< critical resolved stress
crss_back_ID, & !< critical resolved back stress
sense_ID, & !< sense of acting shear stress (-1 or +1)
chi0_ID, & !< backstress at last switch of stress sense (positive?)
gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?)
accshear_ID, &
sumGamma_ID, &
shearrate_ID, &
resolvedstress_ID
end enum
type, private :: tParameters !< container type for internal constitutive parameters
integer(kind(undefined_ID)), dimension(:), allocatable, private :: &
outputID !< ID of each post result output
real(pReal) :: &
gdot0, & !< reference shear strain rate for slip (input parameter)
n_slip, & !< stress exponent for slip (input parameter)
aTolResistance, &
aTolShear
real(pReal), dimension(:), allocatable, private :: &
crss0, & !< initial critical shear stress for slip (input parameter, per family)
theta0, & !< initial hardening rate of forward stress for each slip
theta1, & !< asymptotic hardening rate of forward stress for each slip >
theta0_b, & !< initial hardening rate of back stress for each slip >
theta1_b, & !< asymptotic hardening rate of back stress for each slip >
tau1, &
tau1_b, &
interaction_slipslip, & !< latent hardening matrix
nonSchmidCoeff
real(pReal), dimension(:,:), allocatable, private :: &
hardeningMatrix_SlipSlip
end type
type, private :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress
crss_back, & !< critical resolved back stress
sense, & !< sense of acting shear stress (-1 or +1)
chi0, & !< backstress at last switch of stress sense
gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear
real(pReal), pointer, dimension(:) :: & !< scalars along NipcMyInstance
sumGamma !< accumulated shear across all systems
end type
type(tParameters), dimension(:), allocatable, private :: &
param !< containers of constitutive parameters (len Ninstance)
type(tKinehardeningState), allocatable, dimension(:), private :: &
dotState, &
deltaState, &
state, &
state0
public :: &
plastic_kinehardening_init, &
plastic_kinehardening_LpAndItsTangent, &
plastic_kinehardening_dotState, &
plastic_kinehardening_deltaState, &
plastic_kinehardening_postResults
private :: &
plastic_kinehardening_shearRates
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: &
dEq0
use debug, only: &
debug_level, &
debug_constitutive,&
debug_levelBasic
use math, only: &
math_Mandel3333to66, &
math_Voigt66to3333, &
math_expand
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
PLASTICITY_kinehardening_label, &
PLASTICITY_kinehardening_ID, &
phase_plasticity, &
phase_plasticityInstance, &
phase_Noutput, &
material_phase, &
plasticState, &
MATERIAL_partPhase
use lattice
use numerics,only: &
numerics_integrator
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: &
o, j, k, f, &
output_ID, &
phase, &
instance, &
maxNinstance, &
NipcMyPhase, &
Nchunks_SlipSlip = 0_pInt, Nchunks_SlipFamilies = 0_pInt, &
Nchunks_nonSchmid = 0_pInt, &
offset_slip, index_myFamily, index_otherFamily, &
startIndex, endIndex, &
mySize, nSlip, nSlipFamilies, &
sizeDotState, &
sizeState, &
sizeDeltaState
real(pReal), dimension(:), allocatable :: tempPerSlip
character(len=65536) :: &
tag = '', &
line = '', &
extmsg = ''
character(len=64) :: &
outputtag = ''
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(phase_plasticity == PLASTICITY_KINEHARDENING_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a,1x,i5,/)') '# instances:',maxNinstance
allocate(plastic_kinehardening_sizePostResults(maxNinstance), source=0_pInt)
allocate(plastic_kinehardening_sizePostResult(maxval(phase_Noutput),maxNinstance), &
source=0_pInt)
allocate(plastic_kinehardening_output(maxval(phase_Noutput),maxNinstance))
plastic_kinehardening_output = ''
allocate(plastic_kinehardening_Noutput(maxNinstance), source=0_pInt)
allocate(plastic_kinehardening_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt)
allocate(plastic_kinehardening_totalNslip(maxNinstance), source=0_pInt)
allocate(param(maxNinstance)) ! one container of parameters per instance
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase
phase = phase + 1_pInt ! advance phase section counter
if (phase_plasticity(phase) == PLASTICITY_KINEHARDENING_ID) then
instance = phase_plasticityInstance(phase) ! count instances of my constitutive law
Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) ! maximum number of slip families according to lattice type of current phase
Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase))
Nchunks_nonSchmid = lattice_NnonSchmid(phase)
allocate(param(instance)%outputID(phase_Noutput(phase)), source=undefined_ID) ! allocate space for IDs of every requested output
allocate(param(instance)%crss0 (Nchunks_SlipFamilies), source=0.0_pReal)
allocate(param(instance)%tau1 (Nchunks_SlipFamilies), source=0.0_pReal)
allocate(param(instance)%tau1_b (Nchunks_SlipFamilies), source=0.0_pReal)
allocate(param(instance)%theta0 (Nchunks_SlipFamilies), source=0.0_pReal)
allocate(param(instance)%theta1 (Nchunks_SlipFamilies), source=0.0_pReal)
allocate(param(instance)%theta0_b(Nchunks_SlipFamilies), source=0.0_pReal)
allocate(param(instance)%theta1_b(Nchunks_SlipFamilies), source=0.0_pReal)
allocate(param(instance)%interaction_slipslip(Nchunks_SlipSlip), source=0.0_pReal)
allocate(param(instance)%nonSchmidCoeff(Nchunks_nonSchmid), source=0.0_pReal)
if(allocated(tempPerSlip)) deallocate(tempPerSlip)
allocate(tempPerSlip(Nchunks_SlipFamilies))
endif
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_KINEHARDENING_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
outputtag = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
output_ID = undefined_ID
select case(outputtag)
case ('resistance')
output_ID = crss_ID
case ('backstress')
output_ID = crss_back_ID
case ('sense')
output_ID = sense_ID
case ('chi0')
output_ID = chi0_ID
case ('gamma0')
output_ID = gamma0_ID
case ('accumulatedshear')
output_ID = accshear_ID
case ('totalshear')
output_ID = sumGamma_ID
case ('shearrate')
output_ID = shearrate_ID
case ('resolvedstress')
output_ID = resolvedstress_ID
end select
if (output_ID /= undefined_ID) then
plastic_kinehardening_Noutput(instance) = plastic_kinehardening_Noutput(instance) + 1_pInt
plastic_kinehardening_output(plastic_kinehardening_Noutput(instance),instance) = outputtag
param(instance)%outputID (plastic_kinehardening_Noutput(instance)) = output_ID
endif
!--------------------------------------------------------------------------------------------------
! parameters depending on number of slip families
case ('nslip')
if (chunkPos(1) < Nchunks_SlipFamilies + 1_pInt) &
call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_KINEHARDENING_label//')')
if (chunkPos(1) > Nchunks_SlipFamilies + 1_pInt) &
call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_KINEHARDENING_label//')')
Nchunks_SlipFamilies = chunkPos(1) - 1_pInt ! user specified number of (possibly) active slip families (e.g. 6 0 6 --> 3)
do j = 1_pInt, Nchunks_SlipFamilies
plastic_kinehardening_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
enddo
case ('crss0','tau1','tau1_b','theta0','theta1','theta0_b','theta1_b')
tempPerSlip = 0.0_pReal
do j = 1_pInt, Nchunks_SlipFamilies
if (plastic_kinehardening_Nslip(j,instance) > 0_pInt) &
tempPerSlip(j) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
select case(tag)
case ('crss0')
param(instance)%crss0(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('tau1')
param(instance)%tau1(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('tau1_b')
param(instance)%tau1_b(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('theta0')
param(instance)%theta0(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('theta1')
param(instance)%theta1(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('theta0_b')
param(instance)%theta0_b(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies)
case ('theta1_b')
param(instance)%theta1_b(1:Nchunks_SlipFamilies) = tempPerSlip(1:Nchunks_SlipFamilies)
end select
!--------------------------------------------------------------------------------------------------
! parameters depending on number of interactions
case ('interaction_slipslip')
if (chunkPos(1) < 1_pInt + Nchunks_SlipSlip) &
call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_KINEHARDENING_label//')')
do j = 1_pInt, Nchunks_SlipSlip
param(instance)%interaction_slipslip(j) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
case ('nonschmidcoeff')
if (chunkPos(1) < 1_pInt + Nchunks_nonSchmid) &
call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_KINEHARDENING_label//')')
do j = 1_pInt,Nchunks_nonSchmid
param(instance)%nonSchmidCoeff(j) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
!--------------------------------------------------------------------------------------------------
case ('gdot0')
param(instance)%gdot0 = IO_floatValue(line,chunkPos,2_pInt)
case ('n_slip')
param(instance)%n_slip = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_resistance')
param(instance)%aTolResistance = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_shear')
param(instance)%aTolShear = IO_floatValue(line,chunkPos,2_pInt)
case default
end select
endif; endif
enddo parsingFile
!--------------------------------------------------------------------------------------------------
! allocation of variables whose size depends on the total number of active slip systems
allocate(state(maxNinstance))
allocate(state0(maxNinstance))
allocate(dotState(maxNinstance))
allocate(deltaState(maxNinstance))
initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop through all phases in material.config
myPhase2: if (phase_plasticity(phase) == PLASTICITY_KINEHARDENING_ID) then ! only consider my phase
NipcMyPhase = count(material_phase == phase) ! number of IPCs containing my phase
instance = phase_plasticityInstance(phase) ! which instance of my phase
plastic_kinehardening_Nslip(1:lattice_maxNslipFamily,instance) = &
min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active slip systems per family to min of available and requested
plastic_kinehardening_Nslip(1:lattice_maxNslipFamily,instance))
plastic_kinehardening_totalNslip(instance) = sum(plastic_kinehardening_Nslip(:,instance)) ! how many slip systems altogether
nSlipFamilies = count(plastic_kinehardening_Nslip(:,instance) > 0_pInt)
nSlip = plastic_kinehardening_totalNslip(instance) ! total number of active slip systems
!--------------------------------------------------------------------------------------------------
! sanity checks
if (any(plastic_kinehardening_Nslip(1:nSlipFamilies,instance) > 0_pInt &
.and. param(instance)%crss0(1:nSlipFamilies) < 0.0_pReal)) extmsg = trim(extmsg)//' crss0'
if (any(plastic_kinehardening_Nslip(1:nSlipFamilies,instance) > 0_pInt &
.and. param(instance)%tau1(1:nSlipFamilies) <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1'
if (any(plastic_kinehardening_Nslip(1:nSlipFamilies,instance) > 0_pInt &
.and. param(instance)%tau1_b(1:nSlipFamilies) < 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b'
if (param(instance)%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
if (param(instance)%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip'
if (param(instance)%aTolResistance <= 0.0_pReal) param(instance)%aTolResistance = 1.0_pReal ! default absolute tolerance 1 Pa
if (param(instance)%aTolShear <= 0.0_pReal) param(instance)%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6
if (extmsg /= '') then
extmsg = trim(extmsg)//' ('//PLASTICITY_KINEHARDENING_label//')' ! prepare error message identifier
call IO_error(211_pInt,ip=instance,ext_msg=extmsg)
endif
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,plastic_kinehardening_Noutput(instance)
select case(param(instance)%outputID(o))
case(crss_ID, & !< critical resolved stress
crss_back_ID, & !< critical resolved back stress
sense_ID, & !< sense of acting shear stress (-1 or +1)
chi0_ID, & !< backstress at last switch of stress sense
gamma0_ID, & !< accumulated shear at last switch of stress sense
accshear_ID, &
shearrate_ID, &
resolvedstress_ID)
mySize = nSlip
case(sumGamma_ID)
mySize = 1_pInt
case default
end select
outputFound: if (mySize > 0_pInt) then
plastic_kinehardening_sizePostResult(o,instance) = mySize
plastic_kinehardening_sizePostResults(instance) = plastic_kinehardening_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
!--------------------------------------------------------------------------------------------------
! allocate state arrays
sizeDotState = nSlip & !< crss
+ nSlip & !< crss_back
+ nSlip & !< accumulated (absolute) shear
+ 1_pInt !< sum(gamma)
sizeDeltaState = nSlip & !< sense of acting shear stress (-1 or +1)
+ nSlip & !< backstress at last switch of stress sense
+ nSlip !< accumulated shear at last switch of stress sense
sizeState = sizeDotState + sizeDeltaState
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%offsetDeltaState = sizeDotState
plasticState(phase)%sizePostResults = plastic_kinehardening_sizePostResults(instance)
plasticState(phase)%nSlip = nSlip
allocate(plasticState(phase)%state0 ( sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%partionedState0 ( sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%subState0 ( sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%state ( sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%aTolState (sizeDotState), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase), source=0.0_pReal) ! allocate space for deltaState
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NipcMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal)
offset_slip = plasticState(phase)%nSlip+plasticState(phase)%nTwin+2_pInt
plasticState(phase)%slipRate => &
plasticState(phase)%dotState(offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NipcMyPhase)
plasticState(phase)%accumulatedSlip => &
plasticState(phase)%state(offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NipcMyPhase)
allocate(param(instance)%hardeningMatrix_SlipSlip(nSlip,nSlip), source=0.0_pReal)
do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X
index_myFamily = sum(plastic_kinehardening_Nslip(1:f-1_pInt,instance))
do j = 1_pInt,plastic_kinehardening_Nslip(f,instance) ! loop over (active) systems in my family (slip)
do o = 1_pInt,lattice_maxNslipFamily
index_otherFamily = sum(plastic_kinehardening_Nslip(1:o-1_pInt,instance))
do k = 1_pInt,plastic_kinehardening_Nslip(o,instance) ! loop over (active) systems in other family (slip)
param(instance)%hardeningMatrix_SlipSlip(index_myFamily+j,index_otherFamily+k) = &
param(instance)%interaction_SlipSlip(lattice_interactionSlipSlip( &
sum(lattice_NslipSystem(1:f-1,phase))+j, &
sum(lattice_NslipSystem(1:o-1,phase))+k, &
phase))
enddo; enddo
enddo; enddo
!----------------------------------------------------------------------------------------------
!locally define dotState alias
endindex = 0_pInt
o = endIndex ! offset of dotstate index relative to state index
startIndex = endIndex + 1_pInt
endIndex = endIndex + nSlip
state (instance)%crss => plasticState(phase)%state (startIndex :endIndex ,1:NipcMyPhase)
state0 (instance)%crss => plasticState(phase)%state0 (startIndex :endIndex ,1:NipcMyPhase)
dotState(instance)%crss => plasticState(phase)%dotState (startIndex-o:endIndex-o,1:NipcMyPhase)
state0(instance)%crss = spread(math_expand(param(instance)%crss0,&
plastic_kinehardening_Nslip(:,instance)), &
2, NipcMyPhase)
plasticState(phase)%aTolState(startIndex-o:endIndex-o) = param(instance)%aTolResistance
! .............................................
startIndex = endIndex + 1_pInt
endIndex = endIndex + nSlip
state (instance)%crss_back => plasticState(phase)%state (startIndex :endIndex ,1:NipcMyPhase)
state0 (instance)%crss_back => plasticState(phase)%state0 (startIndex :endIndex ,1:NipcMyPhase)
dotState(instance)%crss_back => plasticState(phase)%dotState (startIndex-o:endIndex-o,1:NipcMyPhase)
state0(instance)%crss_back = 0.0_pReal
plasticState(phase)%aTolState(startIndex-o:endIndex-o) = param(instance)%aTolResistance
! .............................................
startIndex = endIndex + 1_pInt
endIndex = endIndex + nSlip
state (instance)%accshear => plasticState(phase)%state (startIndex :endIndex ,1:NipcMyPhase)
state0 (instance)%accshear => plasticState(phase)%state0 (startIndex :endIndex ,1:NipcMyPhase)
dotState(instance)%accshear => plasticState(phase)%dotState (startIndex-o:endIndex-o,1:NipcMyPhase)
state0(instance)%accshear = 0.0_pReal
plasticState(phase)%aTolState(startIndex-o:endIndex-o) = param(instance)%aTolShear
! .............................................
startIndex = endIndex + 1_pInt
endIndex = endIndex + 1_pInt
state (instance)%sumGamma => plasticState(phase)%state (startIndex ,1:NipcMyPhase)
state0 (instance)%sumGamma => plasticState(phase)%state0 (startIndex ,1:NipcMyPhase)
dotState(instance)%sumGamma => plasticState(phase)%dotState (startIndex-o ,1:NipcMyPhase)
state0(instance)%sumGamma = 0.0_pReal
plasticState(phase)%aTolState(startIndex-o:endIndex-o) = param(instance)%aTolShear
!----------------------------------------------------------------------------------------------
!locally define deltaState alias
o = endIndex
! .............................................
startIndex = endIndex + 1_pInt
endIndex = endIndex + nSlip
state (instance)%sense => plasticState(phase)%state (startIndex :endIndex ,1:NipcMyPhase)
state0 (instance)%sense => plasticState(phase)%state0 (startIndex :endIndex ,1:NipcMyPhase)
deltaState(instance)%sense => plasticState(phase)%deltaState(startIndex-o:endIndex-o,1:NipcMyPhase)
state0(instance)%sense = 0.0_pReal
! .............................................
startIndex = endIndex + 1_pInt
endIndex = endIndex + nSlip
state (instance)%chi0 => plasticState(phase)%state (startIndex :endIndex ,1:NipcMyPhase)
state0 (instance)%chi0 => plasticState(phase)%state0 (startIndex :endIndex ,1:NipcMyPhase)
deltaState(instance)%chi0 => plasticState(phase)%deltaState(startIndex-o:endIndex-o,1:NipcMyPhase)
state0(instance)%chi0 = 0.0_pReal
! .............................................
startIndex = endIndex + 1_pInt
endIndex = endIndex + nSlip
state (instance)%gamma0 => plasticState(phase)%state (startIndex :endIndex ,1:NipcMyPhase)
state0 (instance)%gamma0 => plasticState(phase)%state0 (startIndex :endIndex ,1:NipcMyPhase)
deltaState(instance)%gamma0 => plasticState(phase)%deltaState(startIndex-o:endIndex-o,1:NipcMyPhase)
state0(instance)%gamma0 = 0.0_pReal
endif myPhase2
enddo initializeInstances
end subroutine plastic_kinehardening_init
!--------------------------------------------------------------------------------------------------
!> @brief calculation of shear rates (\dot \gamma)
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_shearRates(gdot_pos,gdot_neg,tau_pos,tau_neg, &
Tstar_v,ph,instance,of)
use lattice, only: &
lattice_NslipSystem, &
lattice_Sslip_v, &
lattice_maxNslipFamily, &
lattice_NnonSchmid
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ph, & !< phase ID
instance, & !< instance of that phase
of !< index of phaseMember
real(pReal), dimension(plastic_kinehardening_totalNslip(instance)), intent(out) :: &
gdot_pos, & !< shear rates from positive line segments
gdot_neg, & !< shear rates from negative line segments
tau_pos, & !< shear stress on positive line segments
tau_neg !< shear stress on negative line segments
integer(pInt) :: &
index_myFamily, &
f,i,j,k
j = 0_pInt
slipFamilies: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
slipSystems: do i = 1_pInt,plastic_kinehardening_Nslip(f,instance)
j = j + 1_pInt
tau_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
tau_neg(j) = tau_pos(j)
nonSchmidSystems: do k = 1,lattice_NnonSchmid(ph)
tau_pos(j) = tau_pos(j) + param(instance)%nonSchmidCoeff(k)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+0,index_myFamily+i,ph))
tau_neg(j) = tau_neg(j) + param(instance)%nonSchmidCoeff(k)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo nonSchmidSystems
enddo slipSystems
enddo slipFamilies
gdot_pos = 0.5_pReal * param(instance)%gdot0 * &
(abs(tau_pos-state(instance)%crss_back(:,of))/ &
state(instance)%crss(:,of))**param(instance)%n_slip &
*sign(1.0_pReal,tau_pos-state(instance)%crss_back(:,of))
gdot_neg = 0.5_pReal * param(instance)%gdot0 * &
(abs(tau_neg-state(instance)%crss_back(:,of))/ &
state(instance)%crss(:,of))**param(instance)%n_slip &
*sign(1.0_pReal,tau_neg-state(instance)%crss_back(:,of))
end subroutine plastic_kinehardening_shearRates
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dTstar99, &
Tstar_v,ipc,ip,el)
use prec, only: &
dNeq0
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g
use math, only: &
math_Plain3333to99, &
math_Mandel6to33, &
math_transpose33
use lattice, only: &
lattice_Sslip, & !< schmid matrix
lattice_Sslip_v, &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_NnonSchmid
use material, only: &
phaseAt, phasememberAt, &
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt) :: &
instance, &
index_myFamily, &
f,i,j,k,l,m,n, &
of, &
ph
real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(phaseAt(ipc,ip,el)))) :: &
gdot_pos,gdot_neg, &
tau_pos,tau_neg
real(pReal) :: &
dgdot_dtau_pos,dgdot_dtau_neg
real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor
real(pReal), dimension(3,3,2) :: &
nonSchmid_tensor
ph = phaseAt(ipc,ip,el) !< figures phase for each material point
of = phasememberAt(ipc,ip,el) !< index of the positions of each constituent of material point, phasememberAt is a function in material that helps figure them out
instance = phase_plasticityInstance(ph)
Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal
dLp_dTstar99 = 0.0_pReal
call plastic_kinehardening_shearRates(gdot_pos,gdot_neg,tau_pos,tau_neg, &
Tstar_v,ph,instance,of)
j = 0_pInt ! reading and marking the starting index for each slip family
slipFamilies: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
slipSystems: do i = 1_pInt,plastic_kinehardening_Nslip(f,instance)
j = j + 1_pInt
! build nonSchmid tensor
nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1)
do k = 1,lattice_NnonSchmid(ph)
nonSchmid_tensor(1:3,1:3,1) = &
nonSchmid_tensor(1:3,1:3,1) + param(instance)%nonSchmidCoeff(k) * &
lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph)
nonSchmid_tensor(1:3,1:3,2) = &
nonSchmid_tensor(1:3,1:3,2) + param(instance)%nonSchmidCoeff(k) * &
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph)
enddo
Lp = Lp + (gdot_pos(j)+gdot_neg(j))*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) ! sum of all gdot*SchmidTensor gives Lp
! Calculation of the tangent of Lp ! sensitivity of Lp
if (dNeq0(gdot_pos(j))) then
dgdot_dtau_pos = gdot_pos(j)*param(instance)%n_slip/(tau_pos(j)-state(instance)%crss_back(j,of))
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = &
dLp_dTstar3333(k,l,m,n) + dgdot_dtau_pos*lattice_Sslip(k,l,1,index_myFamily+i,ph)* &
nonSchmid_tensor(m,n,1)
endif
if (dNeq0(gdot_neg(j))) then
dgdot_dtau_neg = gdot_neg(j)*param(instance)%n_slip/(tau_neg(j)-state(instance)%crss_back(j,of))
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = &
dLp_dTstar3333(k,l,m,n) + dgdot_dtau_neg*lattice_Sslip(k,l,1,index_myFamily+i,ph)* &
nonSchmid_tensor(m,n,2)
endif
enddo slipSystems
enddo slipFamilies
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
end subroutine plastic_kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el)
use prec, only: &
dNeq, &
dEq0
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g
use material, only: &
phaseAt, &
phasememberAt, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(phaseAt(ipc,ip,el)))) :: &
gdot_pos,gdot_neg, &
tau_pos,tau_neg, &
sense
integer(pInt) :: &
ph, &
instance, & !< instance of my instance (unique number of my constitutive model)
of, &
j !< shortcut notation for offset position in state array
ph = phaseAt(ipc,ip,el)
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(ph)
call plastic_kinehardening_shearRates(gdot_pos,gdot_neg,tau_pos,tau_neg, &
Tstar_v,ph,instance,of)
sense = merge(state(instance)%sense(:,of), & ! keep existing...
sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined
dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,'(a)') '======= kinehardening delta state ======='
endif
#endif
!--------------------------------------------------------------------------------------------------
! switch in sense of shear?
do j = 1,plastic_kinehardening_totalNslip(instance)
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,'(i2,1x,f7.4,1x,f7.4)') j,sense(j),state(instance)%sense(j,of)
endif
#endif
if (dNeq(sense(j),state(instance)%sense(j,of),0.1_pReal)) then
deltaState(instance)%sense (j,of) = sense(j) - state(instance)%sense(j,of) ! switch sense
deltaState(instance)%chi0 (j,of) = abs(state(instance)%crss_back(j,of)) - state(instance)%chi0(j,of) ! remember current backstress magnitude
deltaState(instance)%gamma0(j,of) = state(instance)%accshear(j,of) - state(instance)%gamma0(j,of) ! remember current accumulated shear
else
deltaState(instance)%sense (j,of) = 0.0_pReal ! no change
deltaState(instance)%chi0 (j,of) = 0.0_pReal
deltaState(instance)%gamma0(j,of) = 0.0_pReal
endif
enddo
end subroutine plastic_kinehardening_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_dotState(Tstar_v,ipc,ip,el)
use lattice, only: &
lattice_Sslip_v, &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_NnonSchmid
use material, only: &
material_phase, &
phaseAt, phasememberAt, &
plasticState, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation, vector form
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element !< microstructure state
integer(pInt) :: &
instance,ph, &
f,i,j,k, &
index_myFamily,index_otherFamily, &
nSlip, &
offset_accshear, &
of
real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_pos,gdot_neg, &
tau_pos,tau_neg
of = phasememberAt(ipc,ip,el)
ph = phaseAt(ipc,ip,el)
instance = phase_plasticityInstance(ph)
nSlip = plastic_kinehardening_totalNslip(instance)
dotState(instance)%sumGamma(of) = 0.0_pReal
call plastic_kinehardening_shearRates(gdot_pos,gdot_neg,tau_pos,tau_neg, &
Tstar_v,ph,instance,of)
j = 0_pInt
slipFamilies: do f = 1_pInt,lattice_maxNslipFamily
slipSystems: do i = 1_pInt,plastic_kinehardening_Nslip(f,instance)
j = j+1_pInt
dotState(instance)%crss(j,of) = & ! evolution of slip resistance j
dot_product(param(instance)%hardeningMatrix_SlipSlip(j,1:nSlip),abs(gdot_pos+gdot_neg)) * &
( param(instance)%theta1(f) + &
(param(instance)%theta0(f) - param(instance)%theta1(f) &
+ param(instance)%theta0(f)*param(instance)%theta1(f)*state(instance)%sumGamma(of)/param(instance)%tau1(f)) &
*exp(-state(instance)%sumGamma(of)*param(instance)%theta0(f)/param(instance)%tau1(f)) & ! V term depending on the harding law
)
dotState(instance)%crss_back(j,of) = & ! evolution of back stress resistance j
state(instance)%sense(j,of)*abs(gdot_pos(j)+gdot_neg(j)) * &
( param(instance)%theta1_b(f) + &
(param(instance)%theta0_b(f) - param(instance)%theta1_b(f) &
+ param(instance)%theta0_b(f)*param(instance)%theta1_b(f)/(param(instance)%tau1_b(f)+state(instance)%chi0(j,of)) &
*(state(instance)%accshear(j,of)-state(instance)%gamma0(j,of))) &
*exp(-(state(instance)%accshear(j,of)-state(instance)%gamma0(j,of)) &
*param(instance)%theta0_b(f)/(param(instance)%tau1_b(f)+state(instance)%chi0(j,of))) &
) ! V term depending on the harding law for back stress
dotState(instance)%accshear(j,of) = abs(gdot_pos(j)+gdot_neg(j))
dotState(instance)%sumGamma(of) = dotState(instance)%sumGamma(of) + dotState(instance)%accshear(j,of)
enddo slipSystems
enddo slipFamilies
end subroutine plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_kinehardening_postResults(Tstar_v,ipc,ip,el)
use material, only: &
material_phase, &
plasticState, &
phaseAt, phasememberAt, &
phase_plasticityInstance
use lattice, only: &
lattice_Sslip_v, &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_NnonSchmid
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element !< microstructure state
real(pReal), dimension(plastic_kinehardening_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
plastic_kinehardening_postResults
integer(pInt) :: &
instance,ph, of, &
nSlip,&
o,f,i,c,j,k, &
index_myFamily
real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_pos,gdot_neg, &
tau_pos,tau_neg
of = phasememberAt(ipc,ip,el)
ph = phaseAt(ipc,ip,el)
instance = phase_plasticityInstance(ph)
nSlip = plastic_kinehardening_totalNslip(instance)
plastic_kinehardening_postResults = 0.0_pReal
c = 0_pInt
call plastic_kinehardening_shearRates(gdot_pos,gdot_neg,tau_pos,tau_neg, &
Tstar_v,ph,instance,of)
outputsLoop: do o = 1_pInt,plastic_kinehardening_Noutput(instance)
select case(param(instance)%outputID(o))
case (crss_ID)
plastic_kinehardening_postResults(c+1_pInt:c+nSlip) = state(instance)%crss(:,of)
c = c + nSlip
case(crss_back_ID)
plastic_kinehardening_postResults(c+1_pInt:c+nSlip) = state(instance)%crss_back(:,of)
c = c + nSlip
case (sense_ID)
plastic_kinehardening_postResults(c+1_pInt:c+nSlip) = state(instance)%sense(:,of)
c = c + nSlip
case (chi0_ID)
plastic_kinehardening_postResults(c+1_pInt:c+nSlip) = state(instance)%chi0(:,of)
c = c + nSlip
case (gamma0_ID)
plastic_kinehardening_postResults(c+1_pInt:c+nSlip) = state(instance)%gamma0(:,of)
c = c + nSlip
case (accshear_ID)
plastic_kinehardening_postResults(c+1_pInt:c+nSlip) = state(instance)%accshear(:,of)
c = c + nSlip
case (sumGamma_ID)
plastic_kinehardening_postResults(c+1_pInt) = state(instance)%sumGamma(of)
c = c + 1_pInt
case (shearrate_ID)
plastic_kinehardening_postResults(c+1_pInt:c+nSlip) = gdot_pos+gdot_neg
c = c + nSlip
case (resolvedstress_ID)
j = 0_pInt
slipFamilies: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
slipSystems: do i = 1_pInt,plastic_kinehardening_Nslip(f,instance)
j = j + 1_pInt
plastic_kinehardening_postResults(c+j) = &
dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
enddo slipSystems
enddo slipFamilies
c = c + nSlip
end select
enddo outputsLoop
end function plastic_kinehardening_postResults
end module plastic_kinehardening

2
src/prec.f90
View File

@ -68,7 +68,7 @@ module prec
nTrans = 0_pInt
logical :: &
nonlocal = .false.
real(pReal), pointer, dimension(:,:), contiguous :: &
real(pReal), pointer, dimension(:,:) :: &
slipRate, & !< slip rate
accumulatedSlip !< accumulated plastic slip
end type

68
src/spectral_damage.f90
View File

@ -4,6 +4,10 @@
!> @brief Spectral solver for nonlocal damage
!--------------------------------------------------------------------------------------------------
module spectral_damage
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
use prec, only: &
pInt, &
pReal
@ -18,7 +22,6 @@ module spectral_damage
implicit none
private
#include <petsc/finclude/petsc.h90>
character (len=*), parameter, public :: &
spectral_damage_label = 'spectraldamage'
@ -46,13 +49,9 @@ module spectral_damage
public :: &
spectral_damage_init, &
spectral_damage_solution, &
spectral_damage_forward, &
spectral_damage_destroy
spectral_damage_forward
external :: &
PETScFinalize, &
MPI_Abort, &
MPI_Bcast, &
MPI_Allreduce
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -79,32 +78,22 @@ subroutine spectral_damage_init()
damage_nonlocal_getMobility
implicit none
integer(pInt), dimension(:), allocatable :: localK
PetscInt, dimension(:), allocatable :: localK
integer(pInt) :: proc
integer(pInt) :: i, j, k, cell
DM :: damage_grid
Vec :: uBound, lBound
PetscErrorCode :: ierr
character(len=100) :: snes_type
external :: &
SNESCreate, &
SNESSetOptionsPrefix, &
DMDACreate3D, &
SNESSetDM, &
DMDAGetCorners, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
SNESSetFromOptions, &
SNESGetType, &
VecSet, &
DMGetGlobalVector, &
DMRestoreGlobalVector, &
SNESVISetVariableBounds
DMDAGetCorners, &
DMDASNESSetFunctionLocal
write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press, '
write(6,'(/,a)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018 '
write(6,'(a,/)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018 '
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
@ -116,21 +105,23 @@ subroutine spectral_damage_init()
do proc = 1, worldsize
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
call DMDACreate3D(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & !< cut off stencil at boundary
DMDA_STENCIL_BOX, & !< Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & !< global grid
1, 1, worldsize, &
1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap)
grid(1),grid(2),localK, & !< local grid
[grid(1)],[grid(2)],localK, & !< local grid
damage_grid,ierr) !< handle, error
CHKERRQ(ierr)
call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da
call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr)
call DMsetUp(damage_grid,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(damage_grid,solution,ierr); CHKERRQ(ierr) !< global solution vector (grid x 1, i.e. every def grad tensor)
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,spectral_damage_formResidual,&
PETSC_NULL_OBJECT,ierr) !< residual vector of same shape as solution vector
PETSC_NULL_SNES,ierr) !< residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional cli arguments
call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional CLI arguments
call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr)
if (trim(snes_type) == 'vinewtonrsls' .or. &
trim(snes_type) == 'vinewtonssls') then
@ -206,8 +197,7 @@ type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadC
external :: &
VecMin, &
VecMax, &
SNESSolve, &
SNESGetConvergedReason
SNESSolve
spectral_damage_solution%converged =.false.
@ -216,7 +206,7 @@ type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call SNESSolve(damage_snes,PETSC_NULL_OBJECT,solution,ierr); CHKERRQ(ierr)
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)
if (reason < 1) then
@ -244,14 +234,12 @@ type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadC
call VecMin(solution,position,minDamage,ierr); CHKERRQ(ierr)
call VecMax(solution,position,maxDamage,ierr); CHKERRQ(ierr)
if (worldrank == 0) then
if (spectral_damage_solution%converged) &
write(6,'(/,a)') ' ... nonlocal damage converged .....................................'
write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',&
minDamage, maxDamage, stagNorm
write(6,'(/,a)') ' ==========================================================================='
flush(6)
endif
end function spectral_damage_solution
@ -362,9 +350,6 @@ subroutine spectral_damage_forward()
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
external :: &
SNESGetDM
if (cutBack) then
damage_current = damage_lastInc
damage_stagInc = damage_lastInc
@ -399,21 +384,4 @@ subroutine spectral_damage_forward()
end subroutine spectral_damage_forward
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
subroutine spectral_damage_destroy()
implicit none
PetscErrorCode :: ierr
external :: &
VecDestroy, &
SNESDestroy
call VecDestroy(solution,ierr); CHKERRQ(ierr)
call SNESDestroy(damage_snes,ierr); CHKERRQ(ierr)
end subroutine spectral_damage_destroy
end module spectral_damage

24
src/spectral_interface.f90
View File

@ -11,9 +11,9 @@
module DAMASK_interface
use prec, only: &
pInt
implicit none
private
#include <petsc/finclude/petscsys.h>
logical, public, protected :: appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding)
integer(pInt), public, protected :: spectralRestartInc = 0_pInt !< Increment at which calculation starts
character(len=1024), public, protected :: &
@ -44,7 +44,13 @@ contains
subroutine DAMASK_interface_init()
use, intrinsic :: &
iso_fortran_env
#include <petsc/finclude/petscsys.h>
#if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR!=9
===================================================================================================
========================= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =========================
===================================================================================================
#endif
use PETScSys
use system_routines, only: &
getHostName
@ -72,11 +78,8 @@ subroutine DAMASK_interface_init()
logical :: error
external :: &
quit,&
MPI_Comm_rank,&
MPI_Comm_size,&
PETScInitialize, &
MPI_Init_Thread, &
MPI_abort
PETScErrorF, & ! is called in the CHKERRQ macro
PETScInitialize
open(6, encoding='UTF-8') ! for special characters in output
@ -91,7 +94,7 @@ subroutine DAMASK_interface_init()
call quit(1_pInt)
endif
#endif
call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code
call PETScInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code
CHKERRQ(ierr) ! this is a macro definition, it is case sensitive
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,ierr);CHKERRQ(ierr)
@ -104,10 +107,6 @@ subroutine DAMASK_interface_init()
write(output_unit,'(a)') ' STDERR != 0'
call quit(1_pInt)
endif
if (PETSC_VERSION_MAJOR /= 3 .or. PETSC_VERSION_MINOR /= 7) then
write(6,'(a,2(i1.1,a))') 'PETSc ',PETSC_VERSION_MAJOR,'.',PETSC_VERSION_MINOR,'.x not supported'
call quit(1_pInt)
endif
else mainProcess
close(6) ! disable output for non-master processes (open 6 to rank specific file for debug)
open(6,file='/dev/null',status='replace') ! close(6) alone will leave some temp files in cwd
@ -525,5 +524,4 @@ pure function IIO_stringPos(string)
end function IIO_stringPos
end module

72
src/spectral_mech_Basic.f90
View File

@ -2,9 +2,13 @@
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Basic scheme PETSc solver
!> @brief Basic scheme solver
!--------------------------------------------------------------------------------------------------
module spectral_mech_basic
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
use prec, only: &
pInt, &
pReal
@ -16,7 +20,6 @@ module spectral_mech_basic
implicit none
private
#include <petsc/finclude/petsc.h90>
character (len=*), parameter, public :: &
DAMASK_spectral_SolverBasic_label = 'basic'
@ -64,13 +67,9 @@ module spectral_mech_basic
public :: &
basic_init, &
basic_solution, &
basic_forward, &
basic_destroy
basic_forward
external :: &
PETScFinalize, &
MPI_Abort, &
MPI_Bcast, &
MPI_Allreduce
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -118,25 +117,18 @@ subroutine basic_init
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: F
integer(pInt), dimension(:), allocatable :: localK
PetscInt, dimension(:), allocatable :: localK
integer(pInt) :: proc
character(len=1024) :: rankStr
external :: &
SNESCreate, &
SNESSetOptionsPrefix, &
DMDACreate3D, &
SNESSetDM, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions
DMDASNESSetFunctionLocal
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:31–45, 2015'
write(6,'(/,a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
@ -159,16 +151,18 @@ subroutine basic_init
grid(1),grid(2),grid(3), & ! global grid
1 , 1, worldsize, &
9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
grid(1),grid(2),localK, & ! local grid
[grid(1)],[grid(2)],localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Basic_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector
call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
call DMsetUp(da,ierr); CHKERRQ(ierr)
call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,Basic_formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESSetConvergenceTest(snes,Basic_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
call SNESsetConvergenceTest(snes,Basic_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
CHKERRQ(ierr)
call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
!--------------------------------------------------------------------------------------------------
! init fields
@ -255,8 +249,7 @@ type(tSolutionState) function basic_solution(incInfoIn,timeinc,timeinc_old,stres
SNESConvergedReason :: reason
external :: &
SNESSolve, &
SNESGetConvergedReason
SNESsolve
incInfo = incInfoIn
@ -276,13 +269,13 @@ type(tSolutionState) function basic_solution(incInfoIn,timeinc,timeinc_old,stres
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr)
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
Basic_solution%converged = reason > 0
basic_solution%converged = reason > 0
basic_solution%iterationsNeeded = totalIter
basic_solution%termIll = terminallyIll
terminallyIll = .false.
@ -334,10 +327,6 @@ subroutine Basic_formResidual(in,x_scal,f_scal,dummy,ierr)
real(pReal), dimension(3,3) :: &
deltaF_aim
external :: &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
@ -551,25 +540,4 @@ subroutine Basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,s
end subroutine Basic_forward
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
subroutine Basic_destroy()
use spectral_utilities, only: &
Utilities_destroy
implicit none
PetscErrorCode :: ierr
external :: &
VecDestroy, &
SNESDestroy, &
DMDestroy
call VecDestroy(solution_vec,ierr); CHKERRQ(ierr)
call SNESDestroy(snes,ierr); CHKERRQ(ierr)
call DMDestroy(da,ierr); CHKERRQ(ierr)
end subroutine Basic_destroy
end module spectral_mech_basic

68
src/spectral_mech_Polarisation.f90
View File

@ -5,6 +5,10 @@
!> @brief Polarisation scheme solver
!--------------------------------------------------------------------------------------------------
module spectral_mech_Polarisation
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
use prec, only: &
pInt, &
pReal
@ -16,7 +20,6 @@ module spectral_mech_Polarisation
implicit none
private
#include <petsc/finclude/petsc.h90>
character (len=*), parameter, public :: &
DAMASK_spectral_solverPolarisation_label = 'polarisation'
@ -70,13 +73,9 @@ module spectral_mech_Polarisation
public :: &
Polarisation_init, &
Polarisation_solution, &
Polarisation_forward, &
Polarisation_destroy
Polarisation_forward
external :: &
PETScFinalize, &
MPI_Abort, &
MPI_Bcast, &
MPI_Allreduce
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -128,25 +127,18 @@ subroutine Polarisation_init
FandF_tau, & ! overall pointer to solution data
F, & ! specific (sub)pointer
F_tau ! specific (sub)pointer
integer(pInt), dimension(:), allocatable :: localK
PetscInt, dimension(:), allocatable :: localK
integer(pInt) :: proc
character(len=1024) :: rankStr
external :: &
SNESCreate, &
SNESSetOptionsPrefix, &
DMDACreate3D, &
SNESSetDM, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions
DMDASNESsetFunctionLocal
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:31–45, 2015'
write(6,'(/,a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
@ -171,16 +163,18 @@ subroutine Polarisation_init
grid(1),grid(2),grid(3), & ! global grid
1 , 1, worldsize, &
18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
grid(1),grid(2),localK, & ! local grid
[grid(1)],[grid(2)],localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector
call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
call DMsetUp(da,ierr); CHKERRQ(ierr)
call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESSetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
call SNESsetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
CHKERRQ(ierr)
call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
!--------------------------------------------------------------------------------------------------
! init fields
@ -280,8 +274,7 @@ type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_ol
SNESConvergedReason :: reason
external :: &
SNESSolve, &
SNESGetConvergedReason
SNESSolve
incInfo = incInfoIn
@ -304,7 +297,7 @@ type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_ol
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr)
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! check convergence
@ -375,10 +368,6 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
integer(pInt) :: &
i, j, k, e
external :: &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber
F => x_scal(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE)
F_tau => x_scal(1:3,1:3,2,&
@ -685,25 +674,4 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
end subroutine Polarisation_forward
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
subroutine Polarisation_destroy()
use spectral_utilities, only: &
Utilities_destroy
implicit none
PetscErrorCode :: ierr
external :: &
VecDestroy, &
SNESDestroy, &
DMDestroy
call VecDestroy(solution_vec,ierr); CHKERRQ(ierr)
call SNESDestroy(snes,ierr); CHKERRQ(ierr)
call DMDestroy(da,ierr); CHKERRQ(ierr)
end subroutine Polarisation_destroy
end module spectral_mech_Polarisation

67
src/spectral_thermal.f90
View File

@ -4,6 +4,10 @@
!> @brief Spectral solver for thermal conduction
!--------------------------------------------------------------------------------------------------
module spectral_thermal
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
use prec, only: &
pInt, &
pReal
@ -18,7 +22,6 @@ module spectral_thermal
implicit none
private
#include <petsc/finclude/petsc.h90>
character (len=*), parameter, public :: &
spectral_thermal_label = 'spectralthermal'
@ -46,13 +49,9 @@ module spectral_thermal
public :: &
spectral_thermal_init, &
spectral_thermal_solution, &
spectral_thermal_forward, &
spectral_thermal_destroy
spectral_thermal_forward
external :: &
PETScFinalize, &
MPI_Abort, &
MPI_Bcast, &
MPI_Allreduce
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -92,22 +91,15 @@ subroutine spectral_thermal_init
PetscErrorCode :: ierr
external :: &
SNESCreate, &
SNESSetOptionsPrefix, &
DMDACreate3D, &
SNESSetDM, &
DMDAGetCorners, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
SNESSetFromOptions
SNESsetOptionsPrefix, &
DMDAgetCorners, &
DMDASNESsetFunctionLocal
mainProcess: if (worldrank == 0_pInt) then
write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press,'
write(6,'(/,a)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -117,21 +109,23 @@ subroutine spectral_thermal_init
do proc = 1, worldsize
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
call DMDACreate3D(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid
1, 1, worldsize, &
1, 0, & ! #dof (temperature field), ghost boundary width (domain overlap)
grid (1),grid(2),localK, & ! local grid
thermal_grid,ierr) ! handle, error
1, 0, & !< #dof (thermal phase field), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & !< local grid
thermal_grid,ierr) !< handle, error
CHKERRQ(ierr)
call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da
call DMsetFromOptions(thermal_grid,ierr); CHKERRQ(ierr)
call DMsetUp(thermal_grid,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(thermal_grid,solution ,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,spectral_thermal_formResidual,&
PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector
PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
!--------------------------------------------------------------------------------------------------
! init fields
@ -207,8 +201,7 @@ type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,load
external :: &
VecMin, &
VecMax, &
SNESSolve, &
SNESGetConvergedReason
SNESSolve
spectral_thermal_solution%converged =.false.
@ -217,7 +210,7 @@ type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,load
params%timeinc = timeinc
params%timeincOld = timeinc_old
call SNESSolve(thermal_snes,PETSC_NULL_OBJECT,solution,ierr); CHKERRQ(ierr)
call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr)
if (reason < 1) then
@ -247,14 +240,12 @@ type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,load
call VecMin(solution,position,minTemperature,ierr); CHKERRQ(ierr)
call VecMax(solution,position,maxTemperature,ierr); CHKERRQ(ierr)
if (worldrank == 0) then
if (spectral_thermal_solution%converged) &
write(6,'(/,a)') ' ... thermal conduction converged ..................................'
write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',&
minTemperature, maxTemperature, stagNorm
write(6,'(/,a)') ' ==========================================================================='
flush(6)
endif
end function spectral_thermal_solution
@ -361,9 +352,6 @@ subroutine spectral_thermal_forward()
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
external :: &
SNESGetDM
if (cutBack) then
temperature_current = temperature_lastInc
temperature_stagInc = temperature_lastInc
@ -403,21 +391,4 @@ subroutine spectral_thermal_forward()
end subroutine spectral_thermal_forward
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
subroutine spectral_thermal_destroy()
implicit none
PetscErrorCode :: ierr
external :: &
VecDestroy, &
SNESDestroy
call VecDestroy(solution,ierr); CHKERRQ(ierr)
call SNESDestroy(thermal_snes,ierr); CHKERRQ(ierr)
end subroutine spectral_thermal_destroy
end module spectral_thermal

55
src/spectral_utilities.f90
View File

@ -5,6 +5,8 @@
!--------------------------------------------------------------------------------------------------
module spectral_utilities
use, intrinsic :: iso_c_binding
#include <petsc/finclude/petscsys.h>
use PETScSys
use prec, only: &
pReal, &
pInt
@ -13,7 +15,6 @@ module spectral_utilities
implicit none
private
#include <petsc/finclude/petscsys.h>
include 'fftw3-mpi.f03'
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
@ -139,7 +140,6 @@ module spectral_utilities
utilities_constitutiveResponse, &
utilities_calculateRate, &
utilities_forwardField, &
utilities_destroy, &
utilities_updateIPcoords, &
FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
@ -147,6 +147,8 @@ module spectral_utilities
FIELD_DAMAGE_ID
private :: &
utilities_getFreqDerivative
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -195,12 +197,6 @@ subroutine utilities_init()
geomSize
implicit none
external :: &
PETScOptionsClear, &
PETScOptionsInsertString, &
MPI_Abort
PetscErrorCode :: ierr
integer(pInt) :: i, j, k
integer(pInt), dimension(3) :: k_s
@ -214,10 +210,12 @@ subroutine utilities_init()
scalarSize = 1_C_INTPTR_T, &
vecSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T
external :: &
PetscOptionsInsertString
write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>'
write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity, 46:37–53, 2013'
write(6,'(/,a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2012.09.012'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
@ -232,13 +230,13 @@ subroutine utilities_init()
trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.config '; flush(6)
call PetscOptionsClear(PETSC_NULL_OBJECT,ierr)
call PETScOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(PETSCDEBUG),ierr)
if(debugPETSc) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(petsc_defaultOptions),ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(petsc_options),ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
grid1Red = grid(1)/2_pInt + 1_pInt
@ -632,9 +630,6 @@ real(pReal) function utilities_divergenceRMS()
integer(pInt) :: i, j, k, ierr
complex(pReal), dimension(3) :: rescaledGeom
external :: &
MPI_Allreduce
write(6,'(/,a)') ' ... calculating divergence ................................................'
flush(6)
@ -686,9 +681,6 @@ real(pReal) function utilities_curlRMS()
complex(pReal), dimension(3,3) :: curl_fourier
complex(pReal), dimension(3) :: rescaledGeom
external :: &
MPI_Allreduce
write(6,'(/,a)') ' ... calculating curl ......................................................'
flush(6)
@ -962,9 +954,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
external :: &
MPI_Allreduce
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
flush(6)
@ -1081,9 +1070,6 @@ function utilities_forwardField(timeinc,field_lastInc,rate,aim)
real(pReal), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
PetscErrorCode :: ierr
external :: &
MPI_Allreduce
utilities_forwardField = field_lastInc + rate*timeinc
if (present(aim)) then !< correct to match average
fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt
@ -1175,8 +1161,6 @@ subroutine utilities_updateIPcoords(F)
integer(pInt) :: i, j, k, m, ierr
real(pReal), dimension(3) :: step, offset_coords
real(pReal), dimension(3,3) :: Favg
external &
MPI_Bcast
!--------------------------------------------------------------------------------------------------
! integration in Fourier space
@ -1215,21 +1199,4 @@ subroutine utilities_updateIPcoords(F)
end subroutine utilities_updateIPcoords
!--------------------------------------------------------------------------------------------------
!> @brief cleans up
!--------------------------------------------------------------------------------------------------
subroutine utilities_destroy()
implicit none
call fftw_destroy_plan(planTensorForth)
call fftw_destroy_plan(planTensorBack)
call fftw_destroy_plan(planVectorForth)
call fftw_destroy_plan(planVectorBack)
call fftw_destroy_plan(planScalarForth)
call fftw_destroy_plan(planScalarBack)
end subroutine utilities_destroy
end module spectral_utilities