Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Martin Diehl 2016-07-30 08:20:10 +02:00
commit 707b9ba146
31 changed files with 358 additions and 785 deletions

11
CONFIG
View File

@ -1,10 +1,11 @@
# "set"-syntax needed only for tcsh (but works with bash and zsh) # "set"-syntax needed only for tcsh (but works with bash and zsh)
# DAMASK_ROOT will be expanded # 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 DAMASK_NUM_THREADS = 4
set MARC_VERSION=2015
set ABAQUS_VERSION=6.14-5 set MSC_ROOT = /opt/MSC
set MARC_VERSION = 2015
set ABAQUS_VERSION = 6.14-5

View File

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

View File

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

View File

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

View File

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

View File

@ -518,8 +518,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
nCryst, & nCryst, &
numerics_integrator, & numerics_integrator, &
numerics_integrationMode, & numerics_integrationMode, &
numerics_timeSyncing, & numerics_timeSyncing
analyticJaco
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
@ -582,23 +581,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
invFp, & ! inverse of the plastic deformation gradient invFp, & ! inverse of the plastic deformation gradient
Fe_guess, & ! guess for elastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient
Tstar ! 2nd Piola-Kirchhoff stress tensor 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) :: & integer(pInt) :: &
NiterationCrystallite, & ! number of iterations in crystallite loop NiterationCrystallite, & ! number of iterations in crystallite loop
c, & !< counter in integration point component loop c, & !< counter in integration point component loop
@ -1137,371 +1119,115 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! --+>> STIFFNESS CALCULATION <<+-- ! --+>> STIFFNESS CALCULATION <<+--
computeJacobian: if(updateJaco) then 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,& call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), &
!$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNcomponents,error) crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration
elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
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_LiAndItsTangent(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), & temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e), & crystallite_invFi(1:3,1:3,c,i,e)))
c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration rhs_3333 = 0.0_pReal
if (sum(abs(dLidS)) < tol_math_check) then forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dFidS = 0.0_pReal rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33)
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
call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), & temp_3333 = 0.0_pReal
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS 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), & temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFi(1:3,1:3,c,i,e))) crystallite_invFp(1:3,1:3,c,i,e)), &
rhs_3333 = 0.0_pReal math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & 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_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 lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + &
temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & math_mul3333xx3333(dSdFi,dFidS)
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_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & call math_invert(9_pInt,math_identity2nd(9_pInt)+math_Plain3333to99(lhs_3333),temp_99,error)
crystallite_invFp(1:3,1:3,c,i,e)), & if (error) then
math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))) call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, &
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & ext_msg='inversion error in analytic tangent calculation')
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)) 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) + & dFpinvdF = 0.0_pReal
math_mul3333xx3333(dSdFi,dFidS) 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) crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal
if (error) then temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
call IO_warning(warning_ID=600_pInt,el=e,ip=i,g=c, & math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), &
ext_msg='inversion error in analytic tangent calculation') math_transpose33(crystallite_invFp(1:3,1:3,c,i,e))))
dSdF = rhs_3333 forall(p=1_pInt:3_pInt) &
else crystallite_dPdF(p,1:3,p,1:3,c,i,e) = math_transpose33(temp_33)
dSdF = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333)
endif
dFpinvdF = 0.0_pReal temp_33 = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), &
temp_3333 = math_mul3333xx3333(dLpdS,dSdF) math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e)* & 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_inv33(crystallite_subFp0(1:3,1:3,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)
math_mul33x33(temp_3333(1:3,1:3,p,o), &
crystallite_invFi(1:3,1:3,c,i,e)))
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), & crystallite_invFp(1:3,1:3,c,i,e))
math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), & forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))) crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
forall(p=1_pInt:3_pInt) & math_mul33x33(math_mul33x33(temp_33,dSdF(1:3,1:3,p,o)), &
crystallite_dPdF(p,1:3,p,1:3,c,i,e) = math_transpose33(temp_33) 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)), & temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
math_transpose33(crystallite_invFp(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) & math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)))
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
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,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), & enddo; enddo
crystallite_invFp(1:3,1:3,c,i,e)) enddo elementLooping6
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & !$OMP END PARALLEL DO
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
endif computeJacobian endif computeJacobian
!why not OMP? !why not OMP?

View File

@ -1,126 +0,0 @@
! $Id$
! -*- f90 -*-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Note: the syntax of this file is case sensitive.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This file was auto-generated with f2py (version:2_5972).
! See http://cens.ioc.ee/projects/f2py2e/
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The auto-generated file is quite heavily corrected
! For modifying, notice the following hints:
! - if the dimension of an array depend on a array that is itself an input, use the C-Syntax: (1) becomes [0] etc.
! - be sure that the precision defined is integer, real*8, and complex*16
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
python module core ! in
interface ! in :core
module prec
subroutine prec_init
end subroutine prec_init
end module prec
module damask_interface ! in :damask_interface:DAMASK_spectral_interface.f90
subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn) ! in :damask_interface:DAMASK_spectral_interface.f90
character(len=1024), intent(in) :: loadcaseParameterIn
character(len=1024), intent(in) :: geometryParameterIn
end subroutine DAMASK_interface_init
end module damask_interface
module io
subroutine IO_init
end subroutine IO_init
end module io
module numerics
subroutine numerics_init
end subroutine numerics_init
end module numerics
module debug
subroutine debug_init
end subroutine debug_init
end module debug
module math ! in :math:math.f90
subroutine math_init
end subroutine math_init
function math_tensorAvg(field) ! in :math:math.f90
! input variables
real*8 dimension(:,:,:,:,:), intent(in), :: field
! function definition
real*8 dimension(3,3), :: math_tensorAvg
end function math_tensorAvg
end module math
module fesolving
subroutine FE_init
end subroutine FE_init
end module fesolving
module mesh ! in :mesh:mesh.f90
subroutine mesh_init(ip,element)
integer, parameter :: ip = 1
integer, parameter :: element = 1
end subroutine mesh_init
function mesh_nodesAroundCentres(gDim,Favg,centres) ! in :mesh:mesh.f90
real*8, dimension(:,:,:,:), intent(in) :: centres
real*8, dimension(3), intent(in) :: gDim
real*8, dimension(3,3), intent(in) :: Favg
real*8, dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1), depend(centres) :: mesh_nodesAroundCentres
real*8, dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1), depend(centres) :: wrappedCentres
end function mesh_nodesAroundCentres
function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) ! in :mesh:mesh.f90
real*8, dimension(:,:,:,:,:), intent(in) :: F
real*8, dimension(3), intent(in) :: gDim
real*8, dimension(3,3), intent(in), optional :: FavgIn = -1.0
real*8, dimension(3), intent(in), optional :: scalingIn = -1.0
real*8, dimension(3,size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_deformedCoordsFFT
end function mesh_deformedCoordsFFT
function mesh_volumeMismatch(gDim,F,nodes) ! in :mesh:mesh.f90
real*8, dimension(:,:,:,:,:), intent(in) :: F
real*8, dimension(:,:,:,:), intent(in) :: nodes
real*8, dimension(3), intent(in) :: gDim
real*8, dimension(size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_volumeMismatch
end function mesh_volumeMismatch
function mesh_shapeMismatch(gDim,F,nodes,centres) ! in :mesh:mesh.f90
real*8, dimension(:,:,:,:,:), intent(in) :: F
real*8, dimension(:,:,:,:), intent(in) :: nodes
real*8, dimension(:,:,:,:), intent(in) :: centres
real*8, dimension(3), intent(in) :: gDim
real*8, dimension(size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_shapeMismatch
end function mesh_shapeMismatch
function mesh_init_postprocessing(filepath) ! in :mesh:mesh.f90
character(len=*), intent(in) :: filepath
end function mesh_init_postprocessing
function mesh_build_cellnodes(nodes,Ncellnodes) ! in :mesh:mesh.f90
integer, intent(in) :: Ncellnodes
real*8, dimension(3,:), intent(in) :: nodes
real*8, dimension(3,Ncellnodes), depend(Ncellnodes) :: mesh_build_cellnodes
end function mesh_build_cellnodes
function mesh_get_Ncellnodes() ! in :mesh:mesh.f90
integer :: mesh_get_Ncellnodes
end function mesh_get_Ncellnodes
function mesh_get_unitlength() ! in :mesh:mesh.f90
real*8 :: mesh_get_unitlength
end function mesh_get_unitlength
function mesh_get_nodeAtIP(elemtypeFE,ip) ! in :mesh:mesh.f90
character(len=*), intent(in) :: elemtypeFE
integer, intent(in) :: ip
integer :: mesh_get_nodeAtIP
end function mesh_get_nodeAtIP
end module mesh
end interface
end python module core

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -295,8 +295,6 @@ use material, only: phase_plasticity, &
material_phase material_phase
use lattice use lattice
use numerics,only: & use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator numerics_integrator
@ -332,11 +330,9 @@ integer(pInt) :: phase, &
integer(pInt) :: NofMyPhase integer(pInt) :: NofMyPhase
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
maxNinstances = int(count(phase_plasticity == PLASTICITY_NONLOCAL_ID),pInt) 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 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)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,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 if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -0,0 +1,94 @@
[Ni_nonlocal]
elasticity hooke
plasticity nonlocal
/nonlocal/
(output) rho
(output) rho_sgl_mobile
(output) rho_sgl_immobile
(output) rho_sgl_edge_pos
(output) rho_sgl_edge_neg
(output) rho_sgl_screw_pos
(output) rho_sgl_screw_neg
(output) rho_dip_edge
(output) rho_dip_screw
(output) rho_forest
(output) excess_rho_edge
(output) excess_rho_screw
(output) accumulatedshear
(output) shearrate
(output) resolvedstress
(output) resistance
(output) velocity_edge_pos
(output) rho_dot_gen
(output) rho_dot_sgl2dip_edge
(output) rho_dot_sgl2dip_screw
(output) rho_dot_ann_ath
(output) rho_dot_ann_the_edge
(output) rho_dot_ann_the_screw
(output) rho_dot_edgejogs
(output) rho_dot_flux_edge
(output) rho_dot_flux_screw
(output) slipdirection.x
(output) slipdirection.y
(output) slipdirection.z
(output) slipnormal.x
(output) slipnormal.y
(output) slipnormal.z
(output) fluxdensity_edge_pos.x
(output) fluxdensity_edge_pos.y
(output) fluxdensity_edge_pos.z
(output) fluxdensity_edge_neg.x
(output) fluxdensity_edge_neg.y
(output) fluxdensity_edge_neg.z
(output) fluxdensity_screw_pos.x
(output) fluxdensity_screw_pos.y
(output) fluxdensity_screw_pos.z
(output) fluxdensity_screw_neg.x
(output) fluxdensity_screw_neg.y
(output) fluxdensity_screw_neg.z
lattice_structure fcc
Nslip 12 # number of slip systems per family
c11 246.5e9
c12 147.3e9
c44 124.7e9
burgers 2.48e-10 0 0 0 # Burgers vector in m
rhoSglEdgePos0 6e10 # Initial positive edge single dislocation density in m/m**3
rhoSglEdgeNeg0 6e10 # Initial negative edge single dislocation density in m/m**3
rhoSglScrewPos0 6e10 # Initial positive screw single dislocation density in m/m**3
rhoSglScrewNeg0 6e10 # Initial negative screw single dislocation density in m/m**3
rhoDipEdge0 0 # Initial edge dipole dislocation density in m/m**3
rhoDipScrew0 0 # Initial screw dipole dislocation density in m/m**3
rhoSglScatter 0
minimumDipoleHeightEdge 2.6e-9 # 3.0e-9 # minimum distance for stable edge dipoles in m
minimumDipoleHeightScrew 12.0e-9 # 50e-9 # minimum distance for stable screw dipoles in m
lambda0 45 # 33 # prefactor for mean free path
edgeMultiplication 0.1
randomMultiplication 0
atomicVolume 1.2e-29
selfdiffusionPrefactor 1.9e-4 # Gottstein p.168 # prefactor for self-diffusion coefficient
selfdiffusionEnergy 5.1e-19 # Gottstein p.168 # activation energy self-diffusion
solidSolutionEnergy 1.8e-19 # activation energy of solid solution particles in J
solidSolutionConcentration 5e-7 # 1e-7
solidSolutionSize 1.0
peierlsStressEdge 1e5 # Peierls stress for edges in Pa (per slip family)
peierlsStressScrew 1e5 # Peierls stress for screws in Pa (per slip family)
doublekinkWidth 10 # width of double kinks in multiples of burgers vector length b
viscosity 1e-3 # viscosity for dislocation glide in Pa s
p 1 # exponent for thermal barrier profile
q 1 # exponent for thermal barrier profile
attackFrequency 50e9 # attack frequency in Hz
surfaceTransmissivity 1.0 # transmissivity of free surfaces for dislocation flux
grainBoundaryTransmissivity 0.0
aTol_rho 1e100 # absolute tolerance for dislocation density in m/m**3
aTol_shear 1e10 # absolute tolerance for dislocation density in m/m**3
significantRho 1e8 # dislocation density considered relevant in m/m**3
significantN 1
shortRangeStressCorrection 0
CFLfactor 1.1 # safety factor for CFL flux check (numerical parameter)
r 1
interaction_SlipSlip 0 0 0.625 0.07 0.137 0.122 # Dislocation interaction coefficient
linetension 0.8
edgejog 0.01 # 0.2

View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7 #!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,re,sys import os,re,sys
@ -6,10 +6,15 @@ import math
import numpy as np import numpy as np
from optparse import OptionParser from optparse import OptionParser
import damask import damask
import collections
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
def listify(x):
return (x if isinstance(x, collections.Iterable) else [x])
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# MAIN # MAIN
# -------------------------------------------------------------------- # --------------------------------------------------------------------
@ -55,13 +60,9 @@ for i in xrange(len(options.formulas)):
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: try: table = damask.ASCIItable(name = name,
table = damask.ASCIItable(name = name, buffered = False)
buffered = False) except: continue
output = damask.ASCIItable(name = name,
buffered = False)
except:
continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------- # ------------------------------------------ read header -------------------------------------------
@ -129,53 +130,46 @@ for name in filenames:
while outputAlive and table.data_read(): # read next data line of ASCII table while outputAlive and table.data_read(): # read next data line of ASCII table
specials['_row_'] += 1 # count row specials['_row_'] += 1 # count row
output.data_clear()
# ------------------------------------------ calculate one result to get length of labels --------- # ------------------------------------------ calculate one result to get length of labels ---------
if firstLine: if firstLine:
firstLine = False firstLine = False
labelDim = {} resultDim = {}
for label in [x for x in options.labels]: for label in list(options.labels): # iterate over stable copy
labelDim[label] = np.size(eval(evaluator[label])) resultDim[label] = np.size(eval(evaluator[label])) # get dimension of formula[label]
if labelDim[label] == 0: options.labels.remove(label) if resultDim[label] == 0: options.labels.remove(label) # remove label if invalid result
veterans = list(set(options.labels)&set(table.labels(raw=False)+table.labels(raw=True)) ) # intersection of requested and existing
newbies = list(set(options.labels)-set(table.labels(raw=False)+table.labels(raw=True)) ) # requested but not existing
for veteran in list(veterans):
if resultDim[veteran] != table.label_dimension(veteran):
damask.util.croak('{} is ignored due to inconsistent dimension...'.format(veteran))
veterans.remove(veteran) # ignore culprit
for newby in newbies:
table.labels_append(['{}_{}'.format(i+1,newby) for i in xrange(resultDim[newby])]
if resultDim[newby] > 1 else newby)
# ------------------------------------------ assemble header --------------------------------------- # ------------------------------------------ assemble header ---------------------------------------
output.labels_clear() table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
tabLabels = table.labels() table.head_write()
for label in tabLabels:
dim = labelDim[label] if label in options.labels \
else table.label_dimension(label)
output.labels_append(['{}_{}'.format(i+1,label) for i in xrange(dim)] if dim > 1 else label)
for label in options.labels:
if label in tabLabels: continue
output.labels_append(['{}_{}'.format(i+1,label) for i in xrange(labelDim[label])]
if labelDim[label] > 1
else label)
output.info = table.info
output.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
output.head_write()
# ------------------------------------------ process data ------------------------------------------ # ------------------------------------------ process data ------------------------------------------
for label in output.labels(): if options.condition is None or eval(eval(evaluator_condition)): # condition for veteran replacement fulfilled
oldIndices = table.label_indexrange(label) for veteran in veterans: # evaluate formulae that overwrite
Nold = max(1,len(oldIndices)) # Nold could be zero for new columns table.data[table.label_index(veteran):
Nnew = len(output.label_indexrange(label)) table.label_index(veteran)+table.label_dimension(veteran)] = \
output.data_append(eval(evaluator[label]) if label in options.labels and listify(eval(evaluator[veteran]))
(options.condition is None or eval(eval(evaluator_condition)))
else np.tile([table.data[i] for i in oldIndices] for newby in newbies: # evaluate formulae that append
if label in tabLabels table.data_append(listify(eval(evaluator[newby])))
else np.nan,
np.ceil(float(Nnew)/Nold))[:Nnew]) # spread formula result into given number of columns
outputAlive = output.data_write() # output processed line outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------
table.input_close() # close ASCII tables table.close() # close ASCII table
output.close() # close ASCII tables

View File

@ -68,7 +68,7 @@ for name in filenames:
# ------------------------------------------ assemble header -------------------------------------- # ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) 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() table.head_write()
# ------------------------------------------ process data ------------------------------------------ # ------------------------------------------ process data ------------------------------------------

View File

@ -49,7 +49,7 @@ parser.set_defaults(right = (0.0,0.0,0.0))
(options,filename) = parser.parse_args() (options,filename) = parser.parse_args()
if options.format not in outtypes: 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: if options.N < 2:
parser.error('too few steps (need at least 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]: options.trim[0] >= options.trim[1]:
parser.error('invalid trim range (-1 +1).') parser.error('invalid trim range (-1 +1).')
name = options.format if filename == [] \
name = options.format if filename[0] is None\
else filename[0] 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') else open(os.path.basename(filename[0])+extensions[outtypes.index(options.format)],'w')
colorLeft = damask.Color(options.colormodel.upper(), list(options.left)) colorLeft = damask.Color(options.colormodel.upper(), list(options.left))

View File

@ -11,7 +11,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
def mostFrequent(arr): def mostFrequent(arr):
return np.argmax(np.bincount(arr)) return np.argmax(np.bincount(arr.astype('int')))
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
@ -72,7 +72,13 @@ for name in filenames:
# --- do work ------------------------------------------------------------------------------------ # --- do work ------------------------------------------------------------------------------------
microstructure = ndimage.filters.generic_filter(microstructure,mostFrequent,size=(options.stencil,)*3) microstructure = ndimage.filters.generic_filter(microstructure,mostFrequent,size=(options.stencil,)*3).astype('int_')
newInfo = {'microstructures': microstructure.max()}
# --- report ---------------------------------------------------------------------------------------
if ( newInfo['microstructures'] != info['microstructures']):
damask.util.croak('--> microstructures: %i'%newInfo['microstructures'])
info['microstructures'] == newInfo['microstructures']
# --- write header --------------------------------------------------------------------------------- # --- write header ---------------------------------------------------------------------------------