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' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Aritra Chakraborty 2016-07-29 20:42:38 -04:00
parent 304fdf1ebe 1503cc19fd
commit 6a80c71c0b
27 changed files with 219 additions and 614 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
CONFIG
View File

@ -1,10 +1,11 @@
# "set"-syntax needed only for tcsh (but works with bash and zsh)
# DAMASK_ROOT will be expanded
set DAMASK_BIN=${DAMASK_ROOT}/bin
set DAMASK_NUM_THREADS = 4
set DAMASK_BIN = ${DAMASK_ROOT}/bin
set MSC_ROOT=/opt/MSC
set MARC_VERSION=2015
set DAMASK_NUM_THREADS = 4
set ABAQUS_VERSION=6.14-5
set MSC_ROOT = /opt/MSC
set MARC_VERSION = 2015
set ABAQUS_VERSION = 6.14-5

28
DAMASK_env.csh
View File

@ -23,25 +23,29 @@ endif
# according to http://software.intel.com/en-us/forums/topic/501500
# this seems to make sense for the stack size
set FREE=`which free`
set FREE=''
if ( "x$FREE" != "x" ) then
if ( `which free` != "free: Command not found." ) then
set freeMem=`free -k | grep -E '(Mem|Speicher):' | awk '{print $4;}'`
set heap=`expr $freeMem / 2`
set stack=`expr $freeMem / $DAMASK_NUM_THREADS / 2`
set heap=` expr $freeMem / 2`
# http://superuser.com/questions/220059/what-parameters-has-ulimit
limit stacksize $stack # maximum stack size (kB)
limit datasize $heap # maximum heap size (kB)
limit datasize $heap # maximum heap size (kB)
endif
if ( `limit | grep coredumpsize` != "" ) then
limit coredumpsize 0 # prevent core dumping
endif
if ( `limit | grep memoryuse` != "" ) then
limit memoryuse unlimited # maximum physical memory size
endif
if ( `limit | grep vmemoryuse` != "" ) then
limit vmemoryuse unlimited # maximum virtual memory size
endif
limit coredumpsize 2000 # core file size (512-byte blocks)
limit vmemoryuse unlimited # maximum virtual memory size
limit memoryuse unlimited # maximum physical memory size
# disable output in case of scp
if($?prompt) then
if ( $?prompt ) then
echo ''
echo Düsseldorf Advanced Materials Simulation Kit - DAMASK
echo Max-Planck-Institut für Eisenforschung, Düsseldorf
echo Düsseldorf Advanced Materials Simulation Kit --- DAMASK
echo Max-Planck-Institut für Eisenforschung GmbH, Düsseldorf
echo https://damask.mpie.de
echo
echo Using environment with ...
@ -59,8 +63,8 @@ if($?prompt) then
echo "MSC.Marc/Mentat $MSC_ROOT"
endif
echo
echo `limit stacksize`
echo `limit datasize`
echo `limit stacksize`
endif
setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS

54
DAMASK_env.sh
View File

@ -5,8 +5,8 @@
if [ "$OSTYPE" == "linux-gnu" ] || [ "$OSTYPE" == 'linux' ]; then
DAMASK_ROOT=$(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$(dirname $BASH_SOURCE)")
else
[[ "${BASH_SOURCE::1}" == "/" ]] && BASE="" || BASE="`pwd`/"
STAT=$(stat "`dirname $BASE$BASH_SOURCE`")
[[ "${BASH_SOURCE::1}" == "/" ]] && BASE="" || BASE="$(pwd)/"
STAT=$(stat "$(dirname $BASE$BASH_SOURCE)")
DAMASK_ROOT=${STAT##* }
fi
@ -17,16 +17,16 @@ set() {
source $DAMASK_ROOT/CONFIG
unset -f set
# if DAMASK_BIN is present and not in $PATH, add it
if [[ "x$DAMASK_BIN" != "x" && ! `echo ":$PATH:" | grep $DAMASK_BIN:` ]]; then
# add DAMASK_BIN if present but not yet in $PATH
if [[ "x$DAMASK_BIN" != "x" && ! $(echo ":$PATH:" | grep $DAMASK_BIN:) ]]; then
export PATH=$DAMASK_BIN:$PATH
fi
SOLVER=`which DAMASK_spectral 2>/dev/null`
SOLVER=$(which DAMASK_spectral 2>/dev/null)
if [ "x$SOLVER" == "x" ]; then
SOLVER='Not found!'
fi
PROCESSING=`which postResults 2>/dev/null`
PROCESSING=$(which postResults 2>/dev/null)
if [ "x$PROCESSING" == "x" ]; then
PROCESSING='Not found!'
fi
@ -36,25 +36,22 @@ fi
# according to http://software.intel.com/en-us/forums/topic/501500
# this seems to make sense for the stack size
FREE=`which free 2>/dev/null`
FREE=$(which free 2>/dev/null)
if [ "x$FREE" != "x" ]; then
freeMem=`free -k | grep -E '(Mem|Speicher):' | awk '{print $4;}'`
heap=`expr $freeMem / 2`
stack=`expr $freeMem / $DAMASK_NUM_THREADS / 2`
freeMem=$(free -k | grep -E '(Mem|Speicher):' | awk '{print $4;}')
# http://superuser.com/questions/220059/what-parameters-has-ulimit
ulimit -s $stack 2>/dev/null # maximum stack size (kB)
ulimit -d $heap 2>/dev/null # maximum heap size (kB)
ulimit -s $(expr $freeMem / $DAMASK_NUM_THREADS / 2) 2>/dev/null # maximum stack size (kB)
ulimit -d $(expr $freeMem / 2) 2>/dev/null # maximum heap size (kB)
fi
ulimit -c 2000 2>/dev/null # core file size (512-byte blocks)
ulimit -c 0 2>/dev/null # core file size (512-byte blocks)
ulimit -v unlimited 2>/dev/null # maximum virtual memory size
ulimit -m unlimited 2>/dev/null # maximum physical memory size
# disable output in case of scp
if [ ! -z "$PS1" ]; then
echo
echo Düsseldorf Advanced Materials Simulation Kit - DAMASK
echo Max-Planck-Institut für Eisenforschung, Düsseldorf
echo Düsseldorf Advanced Materials Simulation Kit --- DAMASK
echo Max-Planck-Institut für Eisenforschung GmbH, Düsseldorf
echo https://damask.mpie.de
echo
echo Using environment with ...
@ -64,14 +61,29 @@ if [ ! -z "$PS1" ]; then
echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
if [ "x$PETSC_DIR" != "x" ]; then
echo "PETSc location $PETSC_DIR"
[[ `python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR"` == $PETSC_DIR ]] \
|| echo " ~~> "`python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR"`
[[ $(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR") == $PETSC_DIR ]] \
|| echo " ~~> "$(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR")
fi
[[ "x$PETSC_ARCH" != "x" ]] && echo "PETSc architecture $PETSC_ARCH"
[[ "x$PETSC_ARCH" == "x" ]] \
|| echo "PETSc architecture $PETSC_ARCH"
echo "MSC.Marc/Mentat $MSC_ROOT"
echo
echo -n "heap size/MiB "; echo "`ulimit -d`/1024" | bc
echo -n "stack size/MiB "; echo "`ulimit -s`/1024" | bc
echo -n "heap size "
[[ "$(ulimit -d)" == "unlimited" ]] \
&& echo "unlimited" \
|| echo $(python -c \
"import math; \
size=$(( 1024*$(ulimit -d) )); \
print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \
['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0])")
echo -n "stack size "
[[ "$(ulimit -s)" == "unlimited" ]] \
&& echo "unlimited" \
|| echo $(python -c \
"import math; \
size=$(( 1024*$(ulimit -s) )); \
print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \
['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0])")
fi
export DAMASK_NUM_THREADS

54
DAMASK_env.zsh
View File

@ -1,12 +1,7 @@
# sets up an environment for DAMASK on zsh
# usage: source DAMASK_env.zsh
if [ "$OSTYPE" = "linux-gnu" ] || [ "$OSTYPE" = 'linux' ]; then
DAMASK_ROOT=$(readlink -f "`dirname ${(%):-%N}`")
else
print 'Not done yet'
fi
DAMASK_ROOT=${0:a:h}
# defining set() allows to source the same file for tcsh and zsh, with and without space around =
set() {
@ -15,45 +10,36 @@ set() {
source $DAMASK_ROOT/CONFIG
unset -f set
# if DAMASK_BIN is present and not in $PATH, add it
# add DAMASK_BIN if present but not yet in $PATH
MATCH=`echo ":$PATH:" | grep $DAMASK_BIN:`
if [[ ( "x$DAMASK_BIN" != "x" ) && ( "x$MATCH" = "x" ) ]]; then
export PATH=$DAMASK_BIN:$PATH
fi
SOLVER=`which DAMASK_spectral 2>/dev/null`
if [ "x$SOLVER" = "x" ]; then
export SOLVER='Not found!'
fi
PROCESSING=`which postResults 2>/dev/null`
if [ "x$PROCESSING" = "x" ]; then
export PROCESSING='Not found!'
fi
if [ "x$DAMASK_NUM_THREADS" = "x" ]; then
DAMASK_NUM_THREADS=1
fi
# according to http://software.intel.com/en-us/forums/topic/501500
# this seems to make sense for the stack size
FREE=`which free 2>/dev/null`
if [ "x$FREE" != "x" ]; then
if [ "`which free 2>/dev/null`" != "free not found" ]; then
freeMem=`free -k | grep -E '(Mem|Speicher):' | awk '{print $4;}'`
heap=`expr $freeMem / 2`
stack=`expr $freeMem / 2`
# http://superuser.com/questions/220059/what-parameters-has-ulimit
ulimit -s $stack 2>/dev/null # maximum stack size (kB)
ulimit -d $heap 2>/dev/null # maximum heap size (kB)
ulimit -s `expr $freeMem / $DAMASK_NUM_THREADS / 2` 2>/dev/null # maximum stack size (kB)
ulimit -d `expr $freeMem / 2` 2>/dev/null # maximum heap size (kB)
fi
ulimit -c 2000 2>/dev/null # core file size (512-byte blocks)
ulimit -c 0 2>/dev/null # core file size (512-byte blocks)
ulimit -v unlimited 2>/dev/null # maximum virtual memory size
ulimit -m unlimited 2>/dev/null # maximum physical memory size
# disable output in case of scp
if [ ! -z "$PS1" ]; then
echo
echo Düsseldorf Advanced Materials Simulation Kit - DAMASK
echo Max-Planck-Institut für Eisenforschung, Düsseldorf
echo Düsseldorf Advanced Materials Simulation Kit --- DAMASK
echo Max-Planck-Institut für Eisenforschung GmbH, Düsseldorf
echo https://damask.mpie.de
echo
echo Using environment with ...
@ -63,13 +49,29 @@ if [ ! -z "$PS1" ]; then
echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
if [ "x$PETSC_DIR" != "x" ]; then
echo "PETSc location $PETSC_DIR"
[[ `readlink -f $PETSC_DIR` == $PETSC_DIR ]] || echo " ~~> "`readlink -f $PETSC_DIR`
[[ $(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR") == $PETSC_DIR ]] \
|| echo " ~~> "$(python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" "$PETSC_DIR")
fi
[[ "x$PETSC_ARCH" != "x" ]] && echo "PETSc architecture $PETSC_ARCH"
[[ "x$PETSC_ARCH" == "x" ]] \
|| echo "PETSc architecture $PETSC_ARCH"
echo "MSC.Marc/Mentat $MSC_ROOT"
echo
echo -n "heap size/MiB "; echo "`ulimit -d`/1024" | bc
echo -n "stack size/MiB "; echo "`ulimit -s`/1024" | bc
echo -n "heap size "
[[ "$(ulimit -d)" == "unlimited" ]] \
&& echo "unlimited" \
|| echo $(python -c \
"import math; \
size=$(( 1024*$(ulimit -d) )); \
print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \
['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0])")
echo -n "stack size "
[[ "$(ulimit -s)" == "unlimited" ]] \
&& echo "unlimited" \
|| echo $(python -c \
"import math; \
size=$(( 1024*$(ulimit -s) )); \
print '{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \
['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0])")
fi
export DAMASK_NUM_THREADS

2
VERSION
View File

@ -1 +1 @@
v2.0.1
v2.0.1-5-g920cf2c

468
code/crystallite.f90
View File

@ -518,8 +518,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
nCryst, &
numerics_integrator, &
numerics_integrationMode, &
numerics_timeSyncing, &
analyticJaco
numerics_timeSyncing
use debug, only: &
debug_level, &
debug_crystallite, &
@ -582,23 +581,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
invFp, & ! inverse of the plastic deformation gradient
Fe_guess, & ! guess for elastic deformation gradient
Tstar ! 2nd Piola-Kirchhoff stress tensor
real(pReal), allocatable, dimension(:,:,:,:,:,:,:) :: &
dPdF_perturbation1, &
dPdF_perturbation2
real(pReal), allocatable, dimension(:,:,:,:,:) :: &
F_backup, &
Fp_backup, &
InvFp_backup, &
Fi_backup, &
InvFi_backup, &
Fe_backup, &
Lp_backup, &
Li_backup, &
P_backup
real(pReal), allocatable, dimension(:,:,:,:) :: &
Tstar_v_backup
logical, allocatable, dimension(:,:,:) :: &
convergenceFlag_backup
integer(pInt) :: &
NiterationCrystallite, & ! number of iterations in crystallite loop
c, & !< counter in integration point component loop
@ -1137,371 +1119,115 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! --+>> STIFFNESS CALCULATION <<+--
computeJacobian: if(updateJaco) then
jacobianMethod: if (analyticJaco) then
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,&
!$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNcomponents,error)
elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do c = 1_pInt,myNcomponents
call constitutive_TandItsTangent(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent
! --- ANALYTIC JACOBIAN ---
call constitutive_LiAndItsTangent(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e), &
c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration
if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal
else
temp_33 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))
lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) + &
crystallite_subdt(c,i,e)*math_mul33x33(temp_33,dLidFi(1:3,1:3,o,p))
lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) + &
crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e)
rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) - &
crystallite_subdt(c,i,e)*math_mul33x33(temp_33,dLidS(1:3,1:3,o,p))
enddo; enddo
call math_invert(9_pInt,math_Plain3333to99(lhs_3333),temp_99,error)
if (error) then
call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, &
ext_msg='inversion error in analytic tangent calculation')
dFidS = 0.0_pReal
else
dFidS = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333)
endif
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
endif
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,&
!$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNcomponents,error)
elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do c = 1_pInt,myNcomponents
call constitutive_TandItsTangent(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent
call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
call constitutive_LiAndItsTangent(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e), &
c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration
if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal
else
temp_33 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))
lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) + &
crystallite_subdt(c,i,e)*math_mul33x33(temp_33,dLidFi(1:3,1:3,o,p))
lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) + &
crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e)
rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) - &
crystallite_subdt(c,i,e)*math_mul33x33(temp_33,dLidS(1:3,1:3,o,p))
enddo; enddo
call math_invert(9_pInt,math_Plain3333to99(lhs_3333),temp_99,error)
if (error) then
call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, &
ext_msg='inversion error in analytic tangent calculation')
dFidS = 0.0_pReal
else
dFidS = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333)
endif
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
endif
temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
crystallite_invFi(1:3,1:3,c,i,e)))
rhs_3333 = 0.0_pReal
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33)
call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
temp_3333 = 0.0_pReal
temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33,dLpdS(1:3,1:3,p,o)), &
crystallite_invFi(1:3,1:3,c,i,e))
temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
crystallite_invFi(1:3,1:3,c,i,e)))
rhs_3333 = 0.0_pReal
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33)
temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e)), &
math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(temp_33,dLidS(1:3,1:3,p,o))
temp_3333 = 0.0_pReal
temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33,dLpdS(1:3,1:3,p,o)), &
crystallite_invFi(1:3,1:3,c,i,e))
lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + &
math_mul3333xx3333(dSdFi,dFidS)
temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e)), &
math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(temp_33,dLidS(1:3,1:3,p,o))
call math_invert(9_pInt,math_identity2nd(9_pInt)+math_Plain3333to99(lhs_3333),temp_99,error)
if (error) then
call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, &
ext_msg='inversion error in analytic tangent calculation')
dSdF = rhs_3333
else
dSdF = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333)
endif
lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + &
math_mul3333xx3333(dSdFi,dFidS)
dFpinvdF = 0.0_pReal
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e)* &
math_mul33x33(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), &
math_mul33x33(temp_3333(1:3,1:3,p,o), &
crystallite_invFi(1:3,1:3,c,i,e)))
call math_invert(9_pInt,math_identity2nd(9_pInt)+math_Plain3333to99(lhs_3333),temp_99,error)
if (error) then
call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, &
ext_msg='inversion error in analytic tangent calculation')
dSdF = rhs_3333
else
dSdF = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333)
endif
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal
temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e))))
forall(p=1_pInt:3_pInt) &
crystallite_dPdF(p,1:3,p,1:3,c,i,e) = math_transpose33(temp_33)
dFpinvdF = 0.0_pReal
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e)* &
math_mul33x33(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), &
math_mul33x33(temp_3333(1:3,1:3,p,o), &
crystallite_invFi(1:3,1:3,c,i,e)))
temp_33 = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33)
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal
temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e))))
forall(p=1_pInt:3_pInt) &
crystallite_dPdF(p,1:3,p,1:3,c,i,e) = math_transpose33(temp_33)
temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(math_mul33x33(temp_33,dSdF(1:3,1:3,p,o)), &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))
temp_33 = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33)
temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e)), &
math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(temp_33,math_transpose33(dFpinvdF(1:3,1:3,p,o)))
temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(math_mul33x33(temp_33,dSdF(1:3,1:3,p,o)), &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))
temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e)), &
math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(temp_33,math_transpose33(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo
enddo elementLooping6
!$OMP END PARALLEL DO
else jacobianMethod
! --- STANDARD (PERTURBATION METHOD) FOR JACOBIAN ---
numerics_integrationMode = 2_pInt
! --- BACKUP ---
allocate(dPdF_perturbation1(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(dPdF_perturbation2(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(F_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Fp_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(InvFp_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Fi_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(InvFi_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Fe_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Lp_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Li_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(P_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Tstar_v_backup (6, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(convergenceFlag_backup (homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = .false.)
!$OMP PARALLEL DO PRIVATE(myNcomponents)
elementLooping7: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents
plasticState (phaseAt(c,i,e))%state_backup(:,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%state_backup(:,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e))
enddo
plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e))
enddo
F_backup(1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) ! ... and kinematics
Fp_backup(1:3,1:3,c,i,e) = crystallite_Fp(1:3,1:3,c,i,e)
InvFp_backup(1:3,1:3,c,i,e) = crystallite_invFp(1:3,1:3,c,i,e)
Fi_backup(1:3,1:3,c,i,e) = crystallite_Fi(1:3,1:3,c,i,e)
InvFi_backup(1:3,1:3,c,i,e) = crystallite_invFi(1:3,1:3,c,i,e)
Fe_backup(1:3,1:3,c,i,e) = crystallite_Fe(1:3,1:3,c,i,e)
Lp_backup(1:3,1:3,c,i,e) = crystallite_Lp(1:3,1:3,c,i,e)
Li_backup(1:3,1:3,c,i,e) = crystallite_Li(1:3,1:3,c,i,e)
Tstar_v_backup(1:6,c,i,e) = crystallite_Tstar_v(1:6,c,i,e)
P_backup(1:3,1:3,c,i,e) = crystallite_P(1:3,1:3,c,i,e)
convergenceFlag_backup(c,i,e) = crystallite_converged(c,i,e)
enddo; enddo
enddo elementLooping7
!$END PARALLEL DO
! --- CALCULATE STATE AND STRESS FOR PERTURBATION ---
dPdF_perturbation1 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
dPdF_perturbation2 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
pertubationLoop: do perturbation = 1,2 ! forward and backward perturbation
if (iand(pert_method,perturbation) > 0_pInt) then ! mask for desired direction
myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step
do k = 1,3; do l = 1,3 ! ...alter individual components
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)) &
write(6,'(a,2(1x,i1),1x,a,/)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]'
! --- INITIALIZE UNPERTURBED STATE ---
select case(numerics_integrator(numerics_integrationMode))
case(1_pInt)
!why not OMP? ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents
plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%state_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%state_backup(:,phasememberAt(c,i,e))
enddo
plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e))
enddo
crystallite_Fp(1:3,1:3,c,i,e) = Fp_backup(1:3,1:3,c,i,e)
crystallite_invFp(1:3,1:3,c,i,e) = InvFp_backup(1:3,1:3,c,i,e)
crystallite_Fi(1:3,1:3,c,i,e) = Fi_backup(1:3,1:3,c,i,e)
crystallite_invFi(1:3,1:3,c,i,e) = InvFi_backup(1:3,1:3,c,i,e)
crystallite_Fe(1:3,1:3,c,i,e) = Fe_backup(1:3,1:3,c,i,e)
crystallite_Lp(1:3,1:3,c,i,e) = Lp_backup(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = Li_backup(1:3,1:3,c,i,e)
crystallite_Tstar_v(1:6,c,i,e) = Tstar_v_backup(1:6,c,i,e)
enddo; enddo
enddo
case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step
case(4_pInt,5_pInt)
!why not OMP? ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents
plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%subState0(:,phasememberAt(c,i,e))
enddo
plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e))
enddo
crystallite_Fp(1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e)
crystallite_Fi(1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e)
crystallite_Fe(1:3,1:3,c,i,e) = crystallite_subFe0(1:3,1:3,c,i,e)
crystallite_Lp(1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e)
crystallite_Tstar_v(1:6,c,i,e) = crystallite_subTstar0_v(1:6,c,i,e)
enddo; enddo
enddo
end select
! --- PERTURB EITHER FORWARD OR BACKWARD ---
!why not OMP?
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do c = 1,myNcomponents
crystallite_subF(1:3,1:3,c,i,e) = F_backup(1:3,1:3,c,i,e)
crystallite_subF(k,l,c,i,e) = crystallite_subF(k,l,c,i,e) + myPert
crystallite_todo(c,i,e) = crystallite_requested(c,i,e) &
.and. convergenceFlag_backup(c,i,e)
if (crystallite_todo(c,i,e)) crystallite_converged(c,i,e) = .false. ! start out non-converged
enddo; enddo; enddo
select case(numerics_integrator(numerics_integrationMode))
case(1_pInt)
call crystallite_integrateStateFPI()
case(2_pInt)
call crystallite_integrateStateEuler()
case(3_pInt)
call crystallite_integrateStateAdaptiveEuler()
case(4_pInt)
call crystallite_integrateStateRK4()
case(5_pInt)
call crystallite_integrateStateRKCK45()
end select
!why not OMP?
elementLooping8: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
select case(perturbation)
case(1_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. crystallite_converged(c,i,e)) & ! converged state warrants stiffness update
dPdF_perturbation1(1:3,1:3,k,l,c,i,e) = &
(crystallite_P(1:3,1:3,c,i,e) - P_backup(1:3,1:3,c,i,e)) / myPert ! tangent dP_ij/dFg_kl
case(2_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. crystallite_converged(c,i,e)) & ! converged state warrants stiffness update
dPdF_perturbation2(1:3,1:3,k,l,c,i,e) = &
(crystallite_P(1:3,1:3,c,i,e) - P_backup(1:3,1:3,c,i,e)) / myPert ! tangent dP_ij/dFg_kl
end select
enddo elementLooping8
enddo; enddo ! k,l component perturbation loop
endif
enddo pertubationLoop
! --- STIFFNESS ACCORDING TO PERTURBATION METHOD AND CONVERGENCE ---
elementLooping9: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
select case(pert_method)
case(1_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. convergenceFlag_backup(c,i,e)) & ! perturbation mode 1: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,c,i,e)
case(2_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. convergenceFlag_backup(c,i,e)) & ! perturbation mode 2: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = dPdF_perturbation2(1:3,1:3,1:3,1:3,c,i,e)
case(3_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. convergenceFlag_backup(c,i,e)) & ! perturbation mode 3: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,c,i,e) &
+ dPdF_perturbation2(1:3,1:3,1:3,1:3,c,i,e))
end select
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. .not. convergenceFlag_backup(c,i,e)) & ! for any pertubation mode: if central solution did not converge...
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,c,i,e) ! ...use (elastic) fallback
enddo elementLooping9
! --- RESTORE ---
!why not OMP?
elementLooping10: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents
plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%state_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%state_backup(:,phasememberAt(c,i,e))
enddo
plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e))
enddo
crystallite_subF(1:3,1:3,c,i,e) = F_backup(1:3,1:3,c,i,e)
crystallite_Fp(1:3,1:3,c,i,e) = Fp_backup(1:3,1:3,c,i,e)
crystallite_invFp(1:3,1:3,c,i,e) = InvFp_backup(1:3,1:3,c,i,e)
crystallite_Fi(1:3,1:3,c,i,e) = Fi_backup(1:3,1:3,c,i,e)
crystallite_invFi(1:3,1:3,c,i,e) = InvFi_backup(1:3,1:3,c,i,e)
crystallite_Fe(1:3,1:3,c,i,e) = Fe_backup(1:3,1:3,c,i,e)
crystallite_Lp(1:3,1:3,c,i,e) = Lp_backup(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = Li_backup(1:3,1:3,c,i,e)
crystallite_Tstar_v(1:6,c,i,e) = Tstar_v_backup(1:6,c,i,e)
crystallite_P(1:3,1:3,c,i,e) = P_backup(1:3,1:3,c,i,e)
crystallite_converged(c,i,e) = convergenceFlag_backup(c,i,e)
enddo; enddo
enddo elementLooping10
deallocate(dPdF_perturbation1)
deallocate(dPdF_perturbation2)
deallocate(F_backup )
deallocate(Fp_backup )
deallocate(InvFp_backup )
deallocate(Fi_backup )
deallocate(InvFi_backup )
deallocate(Fe_backup )
deallocate(Lp_backup )
deallocate(Li_backup )
deallocate(P_backup )
deallocate(Tstar_v_backup )
deallocate(convergenceFlag_backup)
endif jacobianMethod
enddo; enddo
enddo elementLooping6
!$OMP END PARALLEL DO
endif computeJacobian
!why not OMP?

4
code/numerics.f90
View File

@ -64,7 +64,6 @@ module numerics
charLength = 1.0_pReal, & !< characteristic length scale for gradient problems
residualStiffness = 1.0e-6_pReal !< non-zero residual damage
logical, protected, public :: &
analyticJaco = .true., & !< use analytic Jacobian or perturbation, Default for Spectral solver .true.:
usePingPong = .true., &
numerics_timeSyncing = .false. !< flag indicating if time synchronization in crystallite is used for nonlocal plasticity
@ -315,8 +314,6 @@ subroutine numerics_init
numerics_integrator(1) = IO_intValue(line,chunkPos,2_pInt)
case ('integratorstiffness')
numerics_integrator(2) = IO_intValue(line,chunkPos,2_pInt)
case ('analyticjaco')
analyticJaco = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('usepingpong')
usepingpong = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('timesyncing')
@ -528,7 +525,6 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8)') ' timeSyncing: ',numerics_timeSyncing
write(6,'(a24,1x,L8)') ' analytic Jacobian: ',analyticJaco
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength

12
code/plastic_disloUCLA.f90
View File

@ -152,8 +152,6 @@ subroutine plastic_disloUCLA_init(fileUnit)
MATERIAL_partPhase
use lattice
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -173,11 +171,9 @@ subroutine plastic_disloUCLA_init(fileUnit)
line = ''
real(pReal), dimension(:), allocatable :: tempPerSlip
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOUCLA_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOUCLA_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_DISLOUCLA_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -498,10 +494,6 @@ subroutine plastic_disloUCLA_init(fileUnit)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/plastic_dislotwin.f90
View File

@ -238,8 +238,6 @@ subroutine plastic_dislotwin_init(fileUnit)
MATERIAL_partPhase
use lattice
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -261,11 +259,9 @@ subroutine plastic_dislotwin_init(fileUnit)
line = ''
real(pReal), dimension(:), allocatable :: tempPerSlip, tempPerTwin, tempPerTrans
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_DISLOTWIN_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -930,10 +926,6 @@ subroutine plastic_dislotwin_init(fileUnit)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/plastic_isotropic.f90
View File

@ -94,8 +94,6 @@ subroutine plastic_isotropic_init(fileUnit)
debug_constitutive, &
debug_levelBasic
use numerics, only: &
analyticJaco, &
worldrank, &
numerics_integrator
use math, only: &
math_Mandel3333to66, &
@ -145,11 +143,9 @@ subroutine plastic_isotropic_init(fileUnit)
outputtag = ''
integer(pInt) :: NipcMyPhase
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_ISOTROPIC_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -316,10 +312,6 @@ subroutine plastic_isotropic_init(fileUnit)
allocate(plasticState(phase)%state ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NipcMyPhase),source=0.0_pReal)
endif
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)

9
code/plastic_none.f90
View File

@ -34,7 +34,6 @@ subroutine plastic_none_init
use IO, only: &
IO_timeStamp
use numerics, only: &
worldrank, &
numerics_integrator
use material, only: &
phase_plasticity, &
@ -53,11 +52,9 @@ subroutine plastic_none_init
sizeDotState, &
sizeDeltaState
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONE_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONE_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_none_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -84,11 +81,9 @@ subroutine plastic_none_init
allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase))
allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase))
allocate(plasticState(phase)%state (sizeState,NofMyPhase))
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase))
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase))
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase))
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase))
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase))
allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase))

12
code/plastic_nonlocal.f90
View File

@ -295,8 +295,6 @@ use material, only: phase_plasticity, &
material_phase
use lattice
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
@ -332,11 +330,9 @@ integer(pInt) :: phase, &
integer(pInt) :: NofMyPhase
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstances = int(count(phase_plasticity == PLASTICITY_NONLOCAL_ID),pInt)
if (maxNinstances == 0) return ! we don't have to do anything if there's no instance for this constitutive law
@ -1310,10 +1306,6 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances),
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/plastic_phenoplus.f90
View File

@ -145,8 +145,6 @@ subroutine plastic_phenoplus_init(fileUnit)
MATERIAL_partPhase
use lattice
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -168,11 +166,9 @@ subroutine plastic_phenoplus_init(fileUnit)
line = ''
real(pReal), dimension(:), allocatable :: tempPerSlip
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPLUS_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPLUS_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_PHENOPLUS_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -589,10 +585,6 @@ subroutine plastic_phenoplus_init(fileUnit)
allocate(plasticState(phase)%state ( sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NipcMyPhase),source=0.0_pReal)
endif
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)

12
code/plastic_phenopowerlaw.f90
View File

@ -157,8 +157,6 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
MATERIAL_partPhase
use lattice
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -181,11 +179,9 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
line = ''
real(pReal), dimension(:), allocatable :: tempPerSlip
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -587,10 +583,6 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
allocate(plasticState(phase)%state ( sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NipcMyPhase),source=0.0_pReal)
endif
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)

12
code/plastic_titanmod.f90
View File

@ -216,8 +216,6 @@ subroutine plastic_titanmod_init(fileUnit)
MATERIAL_partPhase
use lattice
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -241,11 +239,9 @@ subroutine plastic_titanmod_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_TITANMOD_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_TITANMOD_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_TITANMOD_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -859,10 +855,6 @@ subroutine plastic_titanmod_init(fileUnit)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

2
code/prec.f90
View File

@ -67,11 +67,9 @@ module prec
real(pReal), allocatable, dimension(:,:) :: &
partionedState0, &
subState0, &
state_backup, &
deltaState, &
previousDotState, & !< state rate of previous xxxx
previousDotState2, & !< state rate two xxxx ago
dotState_backup, & !< backup of state rate
RK4dotState
real(pReal), allocatable, dimension(:,:,:) :: &
RKCK45dotState

12
code/source_damage_anisoBrittle.f90
View File

@ -92,8 +92,6 @@ subroutine source_damage_anisoBrittle_init(fileUnit)
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
use lattice, only: &
lattice_maxNcleavageFamily, &
@ -111,11 +109,9 @@ subroutine source_damage_anisoBrittle_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoBrittle_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoBrittle_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -268,10 +264,6 @@ subroutine source_damage_anisoBrittle_init(fileUnit)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/source_damage_anisoDuctile.f90
View File

@ -96,8 +96,6 @@ subroutine source_damage_anisoDuctile_init(fileUnit)
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
use lattice, only: &
lattice_maxNslipFamily, &
@ -115,11 +113,9 @@ subroutine source_damage_anisoDuctile_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoDuctile_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoDuctile_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_damage_anisoDuctile_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -270,10 +266,6 @@ subroutine source_damage_anisoDuctile_init(fileUnit)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/source_damage_isoBrittle.f90
View File

@ -82,8 +82,6 @@ subroutine source_damage_isoBrittle_init(fileUnit)
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -97,11 +95,9 @@ subroutine source_damage_isoBrittle_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoBrittle_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoBrittle_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -222,10 +218,6 @@ subroutine source_damage_isoBrittle_init(fileUnit)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/source_damage_isoDuctile.f90
View File

@ -82,8 +82,6 @@ subroutine source_damage_isoDuctile_init(fileUnit)
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -97,11 +95,9 @@ subroutine source_damage_isoDuctile_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoDuctile_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoDuctile_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_damage_isoDuctile_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -222,10 +218,6 @@ subroutine source_damage_isoDuctile_init(fileUnit)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/source_thermal_dissipation.f90
View File

@ -68,8 +68,6 @@ subroutine source_thermal_dissipation_init(fileUnit)
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -83,11 +81,9 @@ subroutine source_thermal_dissipation_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -163,10 +159,6 @@ subroutine source_thermal_dissipation_init(fileUnit)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/source_thermal_externalheat.f90
View File

@ -73,8 +73,6 @@ subroutine source_thermal_externalheat_init(fileUnit)
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -89,11 +87,9 @@ subroutine source_thermal_externalheat_init(fileUnit)
line = ''
real(pReal), allocatable, dimension(:,:) :: temp_time, temp_rate
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_thermal_externalheat_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -189,10 +185,6 @@ subroutine source_thermal_externalheat_init(fileUnit)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/source_vacancy_irradiation.f90
View File

@ -70,8 +70,6 @@ subroutine source_vacancy_irradiation_init(fileUnit)
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -85,11 +83,9 @@ subroutine source_vacancy_irradiation_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_irradiation_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_irradiation_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_vacancy_irradiation_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -169,10 +165,6 @@ subroutine source_vacancy_irradiation_init(fileUnit)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/source_vacancy_phenoplasticity.f90
View File

@ -68,8 +68,6 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit)
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -83,11 +81,9 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_phenoplasticity_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_phenoplasticity_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_vacancy_phenoplasticity_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -163,10 +159,6 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

12
code/source_vacancy_thermalfluc.f90
View File

@ -72,8 +72,6 @@ subroutine source_vacancy_thermalfluc_init(fileUnit)
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
@ -87,11 +85,9 @@ subroutine source_vacancy_thermalfluc_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_thermalfluc_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_thermalfluc_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_vacancy_thermalfluc_ID),pInt)
if (maxNinstance == 0_pInt) return
@ -170,10 +166,6 @@ subroutine source_vacancy_thermalfluc_init(fileUnit)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

2
processing/post/addCauchy.py
View File

@ -68,7 +68,7 @@ for name in filenames:
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append(['%i_Cauchy'%(i+1) for i in xrange(9)]) # extend ASCII header with new labels
table.labels_append(['{}_Cauchy'.format(i+1) for i in xrange(9)]) # extend ASCII header with new labels
table.head_write()
# ------------------------------------------ process data ------------------------------------------

7
processing/post/perceptualUniformColorMap.py
View File

@ -49,7 +49,7 @@ parser.set_defaults(right = (0.0,0.0,0.0))
(options,filename) = parser.parse_args()
if options.format not in outtypes:
parser.error('invalid format: "%s" (can be %s).'%(options.format,', '.join(outtypes)))
parser.error('invalid format: "{}" (choices: {}).'.format(options.format,', '.join(outtypes)))
if options.N < 2:
parser.error('too few steps (need at least 2).')
@ -59,10 +59,9 @@ if options.trim[0] < -1.0 or \
options.trim[0] >= options.trim[1]:
parser.error('invalid trim range (-1 +1).')
name = options.format if filename[0] is None\
name = options.format if filename == [] \
else filename[0]
output = sys.stdout if filename[0] is None\
output = sys.stdout if filename == [] \
else open(os.path.basename(filename[0])+extensions[outtypes.index(options.format)],'w')
colorLeft = damask.Color(options.colormodel.upper(), list(options.left))