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

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

This commit is contained in:
Martin Diehl 2016-09-11 19:14:52 +02:00
parent 9d7cad1fa4 77729f39e3
commit fbacdd8e6d
196 changed files with 1648 additions and 2399 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -4,4 +4,5 @@
*.exe *.exe
*.bak *.bak
*~ *~
bin
PRIVATE PRIVATE

11
CONFIG Normal file
View File

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

19
DAMASK_env.bat
View File

@ -1,19 +0,0 @@
:: sets up an environment for DAMASK on Windows
:: usage: call DAMASK_env.bat
@echo off
set LOCATION=%~dp0
set DAMASK_ROOT=%LOCATION%\DAMASK
set DAMASK_NUM_THREADS=2
chcp 1252
Title Düsseldorf Advanced Materials Simulation Kit - DAMASK, MPIE Düsseldorf
echo.
echo Düsseldorf Advanced Materials Simulation Kit - DAMASK
echo Max-Planck-Institut für Eisenforschung, Düsseldorf
echo http://damask.mpie.de
echo.
echo Preparing environment ...
echo DAMASK_ROOT=%DAMASK_ROOT%
echo DAMASK_NUM_THREADS=%DAMASK_NUM_THREADS%
set DAMASK_BIN=%DAMASK_ROOT%\bin
set PATH=%PATH%;%DAMASK_BIN%
set PYTHONPATH=%PYTHONPATH%;%DAMASK_ROOT%\lib

88
DAMASK_env.csh
View File

@ -1,26 +1,68 @@
# sets up an environment for DAMASK on tcsh # sets up an environment for DAMASK on tcsh
# usage: source DAMASK_env.csh # usage: source DAMASK_env.csh
set MAGIG=($_)
set FILENAME=`readlink -f $called[2]` set CALLED=($_)
set LOCATION = `dirname $FILENAME` set DIRNAME=`dirname $CALLED[2]`
setenv DAMASK_ROOT ${LOCATION} set DAMASK_ROOT=`python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" $DIRNAME`
setenv DAMASK_NUM_THREADS 2
# disable output in case of scp source $DAMASK_ROOT/CONFIG
if($?prompt) then
echo # if DAMASK_BIN is present and not in $PATH, add it
echo Düsseldorf Advanced Materials Simulation Kit - DAMASK if ( $?DAMASK_BIN) then
echo Max-Planck-Institut für Eisenforschung, Düsseldorf set MATCH=`echo :${PATH}: | grep ${DAMASK_BIN}:`
echo http://damask.mpie.de if ( "x$MATCH" == "x" ) then
echo set PATH=${DAMASK_BIN}:${PATH}
echo Preparing environment ... endif
echo "DAMASK_ROOT=$DAMASK_ROOT"
echo "DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
endif endif
ulimit -s unlimited
ulimit -c 0 set SOLVER=`which DAMASK_spectral`
ulimit -v unlimited set PROCESSING=`which postResults`
ulimit -m unlimited if ( "x$DAMASK_NUM_THREADS" == "x" ) then
setenv DAMASK_BIN ${DAMASK_ROOT}/bin set DAMASK_NUM_THREADS=1
setenv PATH ${PATH}:${DAMASK_BIN} endif
setenv PYTHONPATH ${PYTHONPATH}:${DAMASK_ROOT}/lib
setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH} # according to http://software.intel.com/en-us/forums/topic/501500
# this seems to make sense for the stack size
if ( `which free` != "free: Command not found." ) then
set freeMem=`free -k | grep -E '(Mem|Speicher):' | awk '{print $4;}'`
set stack=`expr $freeMem / $DAMASK_NUM_THREADS / 2`
set heap=` expr $freeMem / 2`
# http://superuser.com/questions/220059/what-parameters-has-ulimit
limit stacksize $stack # maximum stack size (kB)
limit datasize $heap # maximum heap size (kB)
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
# disable output in case of scp
if ( $?prompt ) then
echo ''
echo Düsseldorf Advanced Materials Simulation Kit --- DAMASK
echo Max-Planck-Institut für Eisenforschung GmbH, Düsseldorf
echo https://damask.mpie.de
echo
echo Using environment with ...
echo "DAMASK $DAMASK_ROOT"
echo "Spectral Solver $SOLVER"
echo "Post Processing $PROCESSING"
echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
if ( $?PETSC_DIR) then
echo "PETSc location $PETSC_DIR"
endif
if ( $?PETSC_ARCH) then
echo "PETSc architecture $PETSC_ARCH"
endif
if ( $?MSC_ROOT) then
echo "MSC.Marc/Mentat $MSC_ROOT"
endif
echo
echo `limit datasize`
echo `limit stacksize`
endif
setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS
setenv PYTHONPATH $DAMASK_ROOT/lib:$PYTHONPATH

80
DAMASK_env.sh
View File

@ -1,77 +1,99 @@
# sets up an environment for DAMASK on bash # sets up an environment for DAMASK on bash
# usage: source DAMASK_env.sh # usage: source DAMASK_env.sh
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
[[ -f $HOME/.damask/damask.conf ]] && source $HOME/.damask/damask.conf || source /etc/damask.conf # defining set() allows to source the same file for tcsh and bash, with and without space around =
set() {
export $1$2$3
}
source $DAMASK_ROOT/CONFIG
unset -f set
# if DAMASK_BIN is present and not in $PATH, add it # add DAMASK_BIN if present but not yet in $PATH
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
if [ "x$DAMASK_NUM_THREADS" == "x" ]; then
DAMASK_NUM_THREADS=1
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 -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 http://damask.mpie.de echo https://damask.mpie.de
echo echo
echo Using environment with ... echo Using environment with ...
echo "DAMASK $DAMASK_ROOT" echo "DAMASK $DAMASK_ROOT"
[[ "x$SOLVER" != "x" ]] && echo "Spectral Solver $SOLVER" echo "Spectral Solver $SOLVER"
[[ "x$PROCESSING" != "x" ]] && echo "Post Processing $PROCESSING" echo "Post Processing $PROCESSING"
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
export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH
for var in BASE STAT SOLVER PROCESSING FREE; do for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN; do
unset "${var}" unset "${var}"
done done
for var in DAMASK IMKL ACML LAPACK MSC FFTW HDF5; do for var in DAMASK MSC; do
unset "${var}_ROOT" unset "${var}_ROOT"
done done
for var in ABAQUS MARC; do
unset "${var}_VERSION"
done

85
DAMASK_env.zsh
View File

@ -1,74 +1,87 @@
# sets up an environment for DAMASK on bash # sets up an environment for DAMASK on zsh
# usage: source DAMASK_env.sh # usage: source DAMASK_env.zsh
if [ "$OSTYPE" = "linux-gnu" ] || [ "$OSTYPE" = 'linux' ]; then DAMASK_ROOT=${0:a:h}
DAMASK_ROOT=$(readlink -f "`dirname ${(%):-%N}`")
else
print 'Not done yet'
fi
[[ -f $HOME/.damask/damask.conf ]] && source $HOME/.damask/damask.conf || source /etc/damask.conf # defining set() allows to source the same file for tcsh and zsh, with and without space around =
set() {
export $1$2$3
}
source $DAMASK_ROOT/CONFIG
unset -f set
# if DAMASK_BIN is present and not in $PATH, add it # add DAMASK_BIN if present but not yet in $PATH
#if [[ [[ "x$DAMASK_BIN" != "x" ]] && ! `echo ":$PATH:" | grep $DAMASK_BIN:` ]]; then MATCH=`echo ":$PATH:" | grep $DAMASK_BIN:`
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 if [ "x$DAMASK_NUM_THREADS" = "x" ]; then
export PROCESSING='Not found!' 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 -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 http://damask.mpie.de echo https://damask.mpie.de
echo echo
echo Using environment with ... echo "Using environment with ..."
echo "DAMASK $DAMASK_ROOT" echo "DAMASK $DAMASK_ROOT"
[[ "x$SOLVER" != "x" ]] && echo "Spectral Solver $SOLVER" echo "Spectral Solver $SOLVER"
[[ "x$PROCESSING" != "x" ]] && echo "Post Processing $PROCESSING" echo "Post Processing $PROCESSING"
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
export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH
for var in BASE STAT SOLVER PROCESSING FREE; do for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN MATCH; do
unset "${var}" unset "${var}"
done done
for var in DAMASK IMKL ACML LAPACK MSC FFTW HDF5; do for var in DAMASK MSC; do
unset "${var}_ROOT" unset "${var}_ROOT"
done done
for var in ABAQUS MARC; do
unset "${var}_VERSION"
done

10
Makefile
View File

@ -17,16 +17,6 @@ FEM:
marc: marc:
@./installation/mods_MarcMentat/apply_DAMASK_modifications.sh ${MAKEFLAGS} @./installation/mods_MarcMentat/apply_DAMASK_modifications.sh ${MAKEFLAGS}
.PHONY: processing
processing:
@if hash cython 2>/dev/null; then \
cd ./lib/damask; \
ln -s orientation.py corientation.pyx; \
CC=gcc python setup_corientation.py build_ext --inplace; \
rm -rv build; \
rm *.c; \
fi
.PHONY: tidy .PHONY: tidy
tidy: tidy:
@$(MAKE) tidy -C code >/dev/null @$(MAKE) tidy -C code >/dev/null

2
VERSION
View File

@ -1 +1 @@
v2.0.0-347-gbe02300 v2.0.1-103-g59e8109

169
code/DAMASK_spectral.f90
View File

@ -81,7 +81,7 @@ program DAMASK_spectral
use spectral_mech_Polarisation use spectral_mech_Polarisation
use spectral_damage use spectral_damage
use spectral_thermal use spectral_thermal
implicit none implicit none
@ -93,9 +93,9 @@ program DAMASK_spectral
logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors
integer(pInt), parameter :: FILEUNIT = 234_pInt !< file unit, DAMASK IO does not support newunit feature integer(pInt), parameter :: FILEUNIT = 234_pInt !< file unit, DAMASK IO does not support newunit feature
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: & integer(pInt) :: &
N_t = 0_pInt, & !< # of time indicators found in load case file N_t = 0_pInt, & !< # of time indicators found in load case file
N_n = 0_pInt, & !< # of increment specifiers found in load case file N_n = 0_pInt, & !< # of increment specifiers found in load case file
N_def = 0_pInt !< # of rate of deformation specifiers found in load case file N_def = 0_pInt !< # of rate of deformation specifiers found in load case file
character(len=65536) :: & character(len=65536) :: &
@ -105,7 +105,7 @@ program DAMASK_spectral
! loop variables, convergence etc. ! loop variables, convergence etc.
real(pReal), dimension(3,3), parameter :: & real(pReal), dimension(3,3), parameter :: &
ones = 1.0_pReal, & ones = 1.0_pReal, &
zeros = 0.0_pReal zeros = 0.0_pReal
integer(pInt), parameter :: & integer(pInt), parameter :: &
subStepFactor = 2_pInt !< for each substep, divide the last time increment by 2.0 subStepFactor = 2_pInt !< for each substep, divide the last time increment by 2.0
real(pReal) :: & real(pReal) :: &
@ -150,18 +150,17 @@ program DAMASK_spectral
MPI_file_get_position, & MPI_file_get_position, &
MPI_file_write, & MPI_file_write, &
MPI_abort, & MPI_abort, &
MPI_finalize, &
MPI_allreduce, & MPI_allreduce, &
PETScFinalize PETScFinalize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt) call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_spectral 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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize field solver information ! initialize field solver information
nActiveFields = 1 nActiveFields = 1
@ -194,19 +193,19 @@ program DAMASK_spectral
call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase
allocate (loadCases(N_n)) ! array of load cases allocate (loadCases(N_n)) ! array of load cases
loadCases%P%myType='p' loadCases%P%myType='p'
do i = 1, size(loadCases) do i = 1, size(loadCases)
allocate(loadCases(i)%ID(nActiveFields)) allocate(loadCases(i)%ID(nActiveFields))
field = 1 field = 1
loadCases(i)%ID(field) = FIELD_MECH_ID ! mechanical active by default loadCases(i)%ID(field) = FIELD_MECH_ID ! mechanical active by default
if (any(thermal_type == THERMAL_conduction_ID)) then ! thermal field active thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
field = field + 1 field = field + 1
loadCases(i)%ID(field) = FIELD_THERMAL_ID loadCases(i)%ID(field) = FIELD_THERMAL_ID
endif endif thermalActive
if (any(damage_type == DAMAGE_nonlocal_ID)) then ! damage field active damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then
field = field + 1 field = field + 1
loadCases(i)%ID(field) = FIELD_DAMAGE_ID loadCases(i)%ID(field) = FIELD_DAMAGE_ID
endif endif damageActive
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -232,12 +231,10 @@ program DAMASK_spectral
endif endif
do j = 1_pInt, 9_pInt do j = 1_pInt, 9_pInt
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a * temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a *
enddo
do j = 1_pInt,9_pInt
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo enddo
loadCases(currentLoadCase)%deformation%maskLogical = & ! logical mask in 3x3 notation loadCases(currentLoadCase)%deformation%maskLogical = & ! logical mask in 3x3 notation
transpose(reshape(temp_maskVector,[ 3,3])) transpose(reshape(temp_maskVector,[ 3,3]))
loadCases(currentLoadCase)%deformation%maskFloat = & ! float (1.0/0.0) mask in 3x3 notation loadCases(currentLoadCase)%deformation%maskFloat = & ! float (1.0/0.0) mask in 3x3 notation
merge(ones,zeros,loadCases(currentLoadCase)%deformation%maskLogical) merge(ones,zeros,loadCases(currentLoadCase)%deformation%maskLogical)
loadCases(currentLoadCase)%deformation%values = math_plain9to33(temp_valueVector) ! values in 3x3 notation loadCases(currentLoadCase)%deformation%values = math_plain9to33(temp_valueVector) ! values in 3x3 notation
@ -245,8 +242,6 @@ program DAMASK_spectral
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
do j = 1_pInt, 9_pInt do j = 1_pInt, 9_pInt
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk
enddo
do j = 1_pInt,9_pInt
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo enddo
loadCases(currentLoadCase)%P%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) loadCases(currentLoadCase)%P%maskLogical = transpose(reshape(temp_maskVector,[ 3,3]))
@ -261,10 +256,10 @@ program DAMASK_spectral
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt) loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt)
loadCases(currentLoadCase)%logscale = 1_pInt loadCases(currentLoadCase)%logscale = 1_pInt
case('freq','frequency','outputfreq') ! frequency of result writings case('freq','frequency','outputfreq') ! frequency of result writings
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt) loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt)
case('r','restart','restartwrite') ! frequency of writing restart information case('r','restart','restartwrite') ! frequency of writing restart information
loadCases(currentLoadCase)%restartfrequency = & loadCases(currentLoadCase)%restartfrequency = &
max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt)) max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt))
case('guessreset','dropguessing') case('guessreset','dropguessing')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('euler') ! rotation of currentLoadCase given in euler angles case('euler') ! rotation of currentLoadCase given in euler angles
@ -273,10 +268,10 @@ program DAMASK_spectral
k = 1_pInt ! assuming keyword indicating degree/radians present k = 1_pInt ! assuming keyword indicating degree/radians present
select case (IO_lc(IO_stringValue(line,chunkPos,i+1_pInt))) select case (IO_lc(IO_stringValue(line,chunkPos,i+1_pInt)))
case('deg','degree') case('deg','degree')
case('rad','radian') ! don't convert from degree to radian case('rad','radian') ! don't convert from degree to radian
l = 0_pInt l = 0_pInt
case default case default
k = 0_pInt k = 0_pInt
end select end select
do j = 1_pInt, 3_pInt do j = 1_pInt, 3_pInt
temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j)
@ -291,7 +286,7 @@ program DAMASK_spectral
loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector)
end select end select
enddo; enddo enddo; enddo
close(FILEUNIT) close(FILEUNIT)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! consistency checks and output of load case ! consistency checks and output of load case
@ -303,14 +298,14 @@ program DAMASK_spectral
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
write(6,'(2x,a)') 'drop guessing along trajectory' write(6,'(2x,a)') 'drop guessing along trajectory'
if (loadCases(currentLoadCase)%deformation%myType=='l') then if (loadCases(currentLoadCase)%deformation%myType == 'l') then
do j = 1_pInt, 3_pInt do j = 1_pInt, 3_pInt
if (any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & if (any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .true.) .and. &
any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .false.)) & any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .false.)) &
errorID = 832_pInt ! each row should be either fully or not at all defined errorID = 832_pInt ! each row should be either fully or not at all defined
enddo enddo
write(6,'(2x,a)') 'velocity gradient:' write(6,'(2x,a)') 'velocity gradient:'
else if (loadCases(currentLoadCase)%deformation%myType=='f') then else if (loadCases(currentLoadCase)%deformation%myType == 'f') then
write(6,'(2x,a)') 'deformation gradient at end of load case:' write(6,'(2x,a)') 'deformation gradient at end of load case:'
else else
write(6,'(2x,a)') 'deformation gradient rate:' write(6,'(2x,a)') 'deformation gradient rate:'
@ -319,12 +314,12 @@ program DAMASK_spectral
if(loadCases(currentLoadCase)%deformation%maskLogical(i,j)) then if(loadCases(currentLoadCase)%deformation%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%deformation%values(i,j) write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%deformation%values(i,j)
else else
write(6,'(2x,12a)',advance='no') ' * ' write(6,'(2x,12a)',advance='no') ' * '
endif endif
enddo; write(6,'(/)',advance='no') enddo; write(6,'(/)',advance='no')
enddo enddo
if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. & if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. &
loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only
if (any(loadCases(currentLoadCase)%P%maskLogical .and. & if (any(loadCases(currentLoadCase)%P%maskLogical .and. &
transpose(loadCases(currentLoadCase)%P%maskLogical) .and. & transpose(loadCases(currentLoadCase)%P%maskLogical) .and. &
reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
@ -334,12 +329,12 @@ program DAMASK_spectral
if(loadCases(currentLoadCase)%P%maskLogical(i,j)) then if(loadCases(currentLoadCase)%P%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%P%values(i,j)*1e-9_pReal write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%P%values(i,j)*1e-9_pReal
else else
write(6,'(2x,12a)',advance='no') ' * ' write(6,'(2x,12a)',advance='no') ' * '
endif endif
enddo; write(6,'(/)',advance='no') enddo; write(6,'(/)',advance='no')
enddo enddo
if (any(abs(math_mul33x33(loadCases(currentLoadCase)%rotation, & if (any(abs(math_mul33x33(loadCases(currentLoadCase)%rotation, &
math_transpose33(loadCases(currentLoadCase)%rotation))-math_I3) >& math_transpose33(loadCases(currentLoadCase)%rotation))-math_I3) > &
reshape(spread(tol_math_check,1,9),[ 3,3]))& reshape(spread(tol_math_check,1,9),[ 3,3]))&
.or. abs(math_det33(loadCases(currentLoadCase)%rotation)) > & .or. abs(math_det33(loadCases(currentLoadCase)%rotation)) > &
1.0_pReal + tol_math_check) errorID = 846_pInt ! given rotation matrix contains strain 1.0_pReal + tol_math_check) errorID = 846_pInt ! given rotation matrix contains strain
@ -360,7 +355,7 @@ program DAMASK_spectral
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing initialization depending on selected solver ! doing initialization depending on selected solver
call Utilities_init() call Utilities_init()
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(1)%ID(field)) select case (loadCases(1)%ID(field))
@ -372,26 +367,26 @@ program DAMASK_spectral
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) & if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) &
call IO_warning(42_pInt, ext_msg='debug Divergence') call IO_warning(42_pInt, ext_msg='debug Divergence')
call AL_init call AL_init
case (DAMASK_spectral_SolverPolarisation_label) case (DAMASK_spectral_SolverPolarisation_label)
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) & if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) &
call IO_warning(42_pInt, ext_msg='debug Divergence') call IO_warning(42_pInt, ext_msg='debug Divergence')
call Polarisation_init call Polarisation_init
case default case default
call IO_error(error_ID = 891, ext_msg = trim(spectral_solver)) call IO_error(error_ID = 891_pInt, ext_msg = trim(spectral_solver))
end select end select
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
call spectral_thermal_init call spectral_thermal_init
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
call spectral_damage_init() call spectral_damage_init()
end select end select
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write header of output file ! write header of output file
if (worldrank == 0) then if (worldrank == 0) then
@ -410,7 +405,7 @@ program DAMASK_spectral
write(resUnit) 'logscales:', loadCases%logscale write(resUnit) 'logscales:', loadCases%logscale
write(resUnit) 'increments:', loadCases%incs ! one entry per LoadCase write(resUnit) 'increments:', loadCases%incs ! one entry per LoadCase
write(resUnit) 'startingIncrement:', restartInc - 1_pInt ! start with writing out the previous inc write(resUnit) 'startingIncrement:', restartInc - 1_pInt ! start with writing out the previous inc
write(resUnit) 'eoh' write(resUnit) 'eoh'
close(resUnit) ! end of header close(resUnit) ! end of header
open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.sta',form='FORMATTED',status='REPLACE') '.sta',form='FORMATTED',status='REPLACE')
@ -429,29 +424,30 @@ program DAMASK_spectral
allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND) allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND)
outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND) outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND)
call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_allreduce') if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_allreduce')
call MPI_file_open(PETSC_COMM_WORLD, & call MPI_file_open(PETSC_COMM_WORLD, &
trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', & trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', &
MPI_MODE_WRONLY + MPI_MODE_APPEND, & MPI_MODE_WRONLY + MPI_MODE_APPEND, &
MPI_INFO_NULL, & MPI_INFO_NULL, &
resUnit, & resUnit, &
ierr) ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_open') if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_open')
call MPI_file_get_position(resUnit,fileOffset,ierr) ! get offset from header call MPI_file_get_position(resUnit,fileOffset,ierr) ! get offset from header
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_get_position') if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_get_position')
fileOffset = fileOffset + sum(outputSize(1:worldrank)) ! offset of my process in file (header + processes before me) fileOffset = fileOffset + sum(outputSize(1:worldrank)) ! offset of my process in file (header + processes before me)
call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr) call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek') if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_seek')
if (.not. appendToOutFile) then ! if not restarting, write 0th increment if (.not. appendToOutFile) then ! if not restarting, write 0th increment
do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, &
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),& call MPI_file_write(resUnit, &
[(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), & reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults,& [(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), &
(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults, &
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr) MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write') if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
if (worldrank == 0) & if (worldrank == 0) &
@ -460,7 +456,7 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loopping over loadcases ! loopping over loadcases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
time0 = time ! currentLoadCase start time time0 = time ! currentLoadCase start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -474,9 +470,9 @@ program DAMASK_spectral
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
else else
if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st currentLoadCase of logarithmic scale else ! not-1st inc of 1st currentLoadCase of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal)) timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal))
endif endif
@ -490,16 +486,16 @@ program DAMASK_spectral
endif endif
timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step
forwarding: if(totalIncsCounter >= restartInc) then forwarding: if (totalIncsCounter >= restartInc) then
stepFraction = 0_pInt stepFraction = 0_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop over sub incs ! loop over sub incs
subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt) subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt)
time = time + timeinc ! forward time time = time + timeinc ! forward time
stepFraction = stepFraction + 1_pInt stepFraction = stepFraction + 1_pInt
remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! report begin of new increment
if (worldrank == 0) then if (worldrank == 0) then
@ -517,7 +513,7 @@ program DAMASK_spectral
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') & ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel '-',stepFraction, '/', subStepFactor**cutBackLevel
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forward fields ! forward fields
@ -543,18 +539,18 @@ program DAMASK_spectral
F_BC = loadCases(currentLoadCase)%deformation, & F_BC = loadCases(currentLoadCase)%deformation, &
P_BC = loadCases(currentLoadCase)%P, & P_BC = loadCases(currentLoadCase)%P, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
end select end select
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
call spectral_thermal_forward (& call spectral_thermal_forward (&
guess,timeinc,timeIncOld,remainingLoadCaseTime) guess,timeinc,timeIncOld,remainingLoadCaseTime)
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
call spectral_damage_forward (& call spectral_damage_forward (&
guess,timeinc,timeIncOld,remainingLoadCaseTime) guess,timeinc,timeIncOld,remainingLoadCaseTime)
end select end select
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve fields ! solve fields
stagIter = 0_pInt stagIter = 0_pInt
@ -570,43 +566,43 @@ program DAMASK_spectral
P_BC = loadCases(currentLoadCase)%P, & P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, & F_BC = loadCases(currentLoadCase)%deformation, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverAL_label) case (DAMASK_spectral_SolverAL_label)
solres(field) = AL_solution (& solres(field) = AL_solution (&
incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
P_BC = loadCases(currentLoadCase)%P, & P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, & F_BC = loadCases(currentLoadCase)%deformation, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverPolarisation_label) case (DAMASK_spectral_SolverPolarisation_label)
solres(field) = Polarisation_solution (& solres(field) = Polarisation_solution (&
incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
P_BC = loadCases(currentLoadCase)%P, & P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, & F_BC = loadCases(currentLoadCase)%deformation, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
end select end select
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
solres(field) = spectral_thermal_solution (& solres(field) = spectral_thermal_solution (&
guess,timeinc,timeIncOld,remainingLoadCaseTime) guess,timeinc,timeIncOld,remainingLoadCaseTime)
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
solres(field) = spectral_damage_solution (& solres(field) = spectral_damage_solution (&
guess,timeinc,timeIncOld,remainingLoadCaseTime) guess,timeinc,timeIncOld,remainingLoadCaseTime)
end select end select
if(.not. solres(field)%converged) exit ! no solution found if (.not. solres(field)%converged) exit ! no solution found
enddo enddo
stagIter = stagIter + 1_pInt stagIter = stagIter + 1_pInt
stagIterate = stagIter < stagItMax .and. & stagIterate = stagIter < stagItMax .and. &
all(solres(:)%converged) .and. & all(solres(:)%converged) .and. &
.not. all(solres(:)%stagConverged) .not. all(solres(:)%stagConverged)
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check solution ! check solution
cutBack = .False. cutBack = .False.
if(solres(1)%termIll .or. .not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if(solres(1)%termIll .or. .not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back if (cutBackLevel < maxCutBack) then ! do cut back
if (worldrank == 0) write(6,'(/,a)') ' cut back detected' if (worldrank == 0) write(6,'(/,a)') ' cut back detected'
@ -619,8 +615,8 @@ program DAMASK_spectral
call IO_warning(850_pInt) call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
elseif (continueCalculation == 1_pInt) then elseif (continueCalculation == 1_pInt) then
guess = .true. ! accept non converged BVP solution guess = .true. ! accept non converged BVP solution
else ! default behavior, exit if spectral solver does not converge else ! default behavior, exit if spectral solver does not converge
call IO_warning(850_pInt) call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
endif endif
@ -632,8 +628,8 @@ program DAMASK_spectral
write(statUnit,*) totalIncsCounter, time, cutBackLevel, & write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
flush(statUnit) flush(statUnit)
endif endif
endif endif
enddo subIncLooping enddo subIncLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
if(all(solres(:)%converged)) then ! report converged inc if(all(solres(:)%converged)) then ! report converged inc
@ -664,11 +660,11 @@ program DAMASK_spectral
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif endif
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! first call to CPFEM_general will write? mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! first call to CPFEM_general will write?
restartWrite = .true. restartWrite = .true.
lastRestartWritten = inc lastRestartWritten = inc
endif endif
else forwarding else forwarding
time = time + timeinc time = time + timeinc
guess = .true. guess = .true.
@ -700,7 +696,7 @@ program DAMASK_spectral
call AL_destroy() call AL_destroy()
case (DAMASK_spectral_SolverPolarisation_label) case (DAMASK_spectral_SolverPolarisation_label)
call Polarisation_destroy() call Polarisation_destroy()
end select end select
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
call spectral_thermal_destroy() call spectral_thermal_destroy()
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
@ -711,6 +707,13 @@ program DAMASK_spectral
call PETScFinalize(ierr); CHKERRQ(ierr) call PETScFinalize(ierr); CHKERRQ(ierr)
#ifdef _OPENMP
call MPI_finalize(i)
if (i /= 0_pInt) then
call IO_error(error_ID=894, el=i, ext_msg="Finalize()")
endif
#endif
if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;) call quit(0_pInt) ! no complains ;)
@ -721,7 +724,7 @@ end program DAMASK_spectral
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief quit subroutine to mimic behavior of FEM solvers !> @brief quit subroutine to mimic behavior of FEM solvers
!> @details exits the Spectral solver and reports time and duration. Exit code 0 signals !> @details exits the Spectral solver and reports time and duration. Exit code 0 signals
!> everything went fine. Exit code 1 signals an error, message according to IO_error. Exit code !> everything went fine. Exit code 1 signals an error, message according to IO_error. Exit code
!> 2 signals no converged solution and increment of last saved restart information is written to !> 2 signals no converged solution and increment of last saved restart information is written to
!> stderr. Exit code 3 signals no severe problems, but some increments did not converge !> stderr. Exit code 3 signals no severe problems, but some increments did not converge
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -741,7 +744,7 @@ subroutine quit(stop_id)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',& write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',& dateAndTime(6),':',&
dateAndTime(7) dateAndTime(7)
if (stop_id == 0_pInt) stop 0 ! normal termination if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id < 0_pInt) then ! terminally ill, restart might help if (stop_id < 0_pInt) then ! terminally ill, restart might help
write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt) write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)

2
code/FEsolving.f90
View File

@ -60,8 +60,6 @@ subroutine FE_init
IO_warning, & IO_warning, &
IO_timeStamp IO_timeStamp
use DAMASK_interface use DAMASK_interface
use numerics, only: &
worldrank
implicit none implicit none
#if defined(Marc4DAMASK) || defined(Abaqus) #if defined(Marc4DAMASK) || defined(Abaqus)

3
code/IO.f90
View File

@ -129,6 +129,7 @@ recursive function IO_read(fileUnit,reset) result(line)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! normal case ! normal case
if (input == '') return ! regular line if (input == '') return ! regular line
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! recursion case ! recursion case
if (stack >= 10_pInt) call IO_error(104_pInt,ext_msg=input) ! recursion limit reached if (stack >= 10_pInt) call IO_error(104_pInt,ext_msg=input) ! recursion limit reached
@ -141,7 +142,7 @@ recursive function IO_read(fileUnit,reset) result(line)
pathOn(stack) = path(1:scan(path,SEP,.true.))//input ! glue include to current file's dir pathOn(stack) = path(1:scan(path,SEP,.true.))//input ! glue include to current file's dir
endif endif
open(newunit=unitOn(stack),iostat=myStat,file=pathOn(stack)) ! open included file open(newunit=unitOn(stack),iostat=myStat,file=pathOn(stack),action='read') ! open included file
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=pathOn(stack)) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=pathOn(stack))
line = IO_read(fileUnit) line = IO_read(fileUnit)

115
code/Makefile
View File

@ -7,7 +7,7 @@ SHELL = /bin/sh
# OPTIONS = standard (alternative): meaning # OPTIONS = standard (alternative): meaning
#------------------------------------------------------------- #-------------------------------------------------------------
# F90 = ifort (gfortran): compiler type, choose Intel or GNU # F90 = ifort (gfortran): compiler type, choose Intel or GNU
# FCOMPILERNAME = name of the compiler executable (if not the same as the ype), e.g. using mpich-g90 instead of ifort # FCOMPILERNAME = name of the compiler executable (if not the same as the type), e.g. using mpich-g90 instead of ifort
# PORTABLE = TRUE (FALSE): decision, if executable is optimized for the machine on which it was built. # PORTABLE = TRUE (FALSE): decision, if executable is optimized for the machine on which it was built.
# OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE,ULTRA): Optimization mode: O2, O0, O3 + further options for most files, O3 + further options for all files # OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE,ULTRA): Optimization mode: O2, O0, O3 + further options for most files, O3 + further options for all files
# OPENMP = TRUE (FALSE): OpenMP multiprocessor support # OPENMP = TRUE (FALSE): OpenMP multiprocessor support
@ -257,10 +257,10 @@ COMPILE_OPTIONS_gfortran :=-DDAMASKVERSION=\"${DAMASKVERSION}\"\
#-Wunsafe-loop-optimizations: warn if the loop cannot be optimized due to nontrivial assumptions. #-Wunsafe-loop-optimizations: warn if the loop cannot be optimized due to nontrivial assumptions.
#-Wstrict-overflow: #-Wstrict-overflow:
DEBUG_OPTIONS_gfortran :=-g\ DEBUG_OPTIONS_gfortran :=-g \
-fbacktrace\ -fbacktrace \
-fdump-core\ -fdump-core \
-fcheck=all\ -fcheck=all \
-ffpe-trap=invalid,zero,overflow -ffpe-trap=invalid,zero,overflow
################################################################################################### ###################################################################################################
@ -300,37 +300,60 @@ COMPILE =$(OPENMP_FLAG_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$
COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90))
################################################################################################### ###################################################################################################
SOURCE_FILES = \ SOURCE_FILES = \
source_thermal_dissipation.o source_thermal_externalheat.o \ source_thermal_dissipation.o \
source_damage_isoBrittle.o source_damage_isoDuctile.o source_damage_anisoBrittle.o source_damage_anisoDuctile.o \ source_thermal_externalheat.o \
source_vacancy_phenoplasticity.o source_vacancy_irradiation.o source_vacancy_thermalfluc.o source_damage_isoBrittle.o \
source_damage_isoDuctile.o \
source_damage_anisoBrittle.o \
source_damage_anisoDuctile.o \
source_vacancy_phenoplasticity.o \
source_vacancy_irradiation.o \
source_vacancy_thermalfluc.o
KINEMATICS_FILES = \ KINEMATICS_FILES = \
kinematics_cleavage_opening.o kinematics_slipplane_opening.o \ kinematics_cleavage_opening.o \
kinematics_slipplane_opening.o \
kinematics_thermal_expansion.o \ kinematics_thermal_expansion.o \
kinematics_vacancy_strain.o kinematics_hydrogen_strain.o kinematics_vacancy_strain.o \
kinematics_hydrogen_strain.o
PLASTIC_FILES = \ PLASTIC_FILES = \
plastic_dislotwin.o plastic_disloUCLA.o plastic_isotropic.o \ plastic_dislotwin.o \
plastic_phenopowerlaw.o plastic_titanmod.o plastic_nonlocal.o plastic_none.o \ plastic_disloUCLA.o \
plastic_isotropic.o \
plastic_phenopowerlaw.o \
plastic_titanmod.o \
plastic_nonlocal.o \
plastic_none.o \
plastic_phenoplus.o plastic_phenoplus.o
THERMAL_FILES = \ THERMAL_FILES = \
thermal_isothermal.o thermal_adiabatic.o thermal_conduction.o thermal_isothermal.o \
thermal_adiabatic.o \
thermal_conduction.o
DAMAGE_FILES = \ DAMAGE_FILES = \
damage_none.o damage_local.o damage_nonlocal.o damage_none.o \
damage_local.o \
damage_nonlocal.o
VACANCYFLUX_FILES = \ VACANCYFLUX_FILES = \
vacancyflux_isoconc.o vacancyflux_isochempot.o vacancyflux_cahnhilliard.o vacancyflux_isoconc.o \
vacancyflux_isochempot.o \
vacancyflux_cahnhilliard.o
POROSITY_FILES = \ POROSITY_FILES = \
porosity_none.o porosity_phasefield.o porosity_none.o \
porosity_phasefield.o
HYDROGENFLUX_FILES = \ HYDROGENFLUX_FILES = \
hydrogenflux_isoconc.o hydrogenflux_cahnhilliard.o hydrogenflux_isoconc.o \
hydrogenflux_cahnhilliard.o
HOMOGENIZATION_FILES = \ HOMOGENIZATION_FILES = \
homogenization_RGC.o homogenization_isostrain.o homogenization_none.o homogenization_RGC.o \
homogenization_isostrain.o \
homogenization_none.o
##################### #####################
# Spectral Solver # Spectral Solver
@ -351,11 +374,28 @@ DAMASK_spectral.o: INTERFACENAME := spectral_interface.f90
SPECTRAL_SOLVER_FILES = spectral_mech_AL.o spectral_mech_Basic.o spectral_mech_Polarisation.o \ SPECTRAL_SOLVER_FILES = spectral_mech_AL.o spectral_mech_Basic.o spectral_mech_Polarisation.o \
spectral_thermal.o spectral_damage.o spectral_thermal.o spectral_damage.o
SPECTRAL_FILES = C_routines.o system_routines.o prec.o DAMASK_interface.o IO.o numerics.o debug.o math.o \ SPECTRAL_FILES = C_routines.o \
FEsolving.o mesh.o material.o lattice.o \ system_routines.o \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \ prec.o \
DAMASK_interface.o \
IO.o \
numerics.o \
debug.o \
math.o \
FEsolving.o \
mesh.o \
material.o \
lattice.o \
$(SOURCE_FILES) \
$(KINEMATICS_FILES) \
$(PLASTIC_FILES) \
constitutive.o \
crystallite.o \ crystallite.o \
$(THERMAL_FILES) $(DAMAGE_FILES) $(VACANCYFLUX_FILES) $(HYDROGENFLUX_FILES) $(POROSITY_FILES) \ $(THERMAL_FILES) \
$(DAMAGE_FILES) \
$(VACANCYFLUX_FILES) \
$(HYDROGENFLUX_FILES) \
$(POROSITY_FILES) \
$(HOMOGENIZATION_FILES) homogenization.o \ $(HOMOGENIZATION_FILES) homogenization.o \
CPFEM2.o \ CPFEM2.o \
spectral_utilities.o \ spectral_utilities.o \
@ -401,14 +441,31 @@ DAMASK_FEM.exe: INCLUDE_DIRS += -I./
FEM_SOLVER_FILES = FEM_mech.o FEM_thermal.o FEM_damage.o FEM_vacancyflux.o FEM_porosity.o FEM_hydrogenflux.o FEM_SOLVER_FILES = FEM_mech.o FEM_thermal.o FEM_damage.o FEM_vacancyflux.o FEM_porosity.o FEM_hydrogenflux.o
FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o numerics.o debug.o math.o \ FEM_FILES = prec.o \
FEsolving.o mesh.o material.o lattice.o \ DAMASK_interface.o \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \ FEZoo.o \
IO.o \
numerics.o \
debug.o \
math.o \
FEsolving.o \
mesh.o \
material.o \
lattice.o \
$(SOURCE_FILES) \
$(KINEMATICS_FILES) \
$(PLASTIC_FILES) \
constitutive.o \
crystallite.o \ crystallite.o \
$(THERMAL_FILES) $(DAMAGE_FILES) $(VACANCYFLUX_FILES) $(HYDROGENFLUX_FILES) $(POROSITY_FILES) \ $(THERMAL_FILES) \
$(DAMAGE_FILES) \
$(VACANCYFLUX_FILES) \
$(HYDROGENFLUX_FILES) \
$(POROSITY_FILES) \
$(HOMOGENIZATION_FILES) homogenization.o \ $(HOMOGENIZATION_FILES) homogenization.o \
CPFEM.o \ CPFEM.o \
FEM_utilities.o $(FEM_SOLVER_FILES) FEM_utilities.o \
$(FEM_SOLVER_FILES)
DAMASK_FEM.exe: DAMASK_FEM_driver.o DAMASK_FEM.exe: DAMASK_FEM_driver.o
$(PREFIX) $(LINKERNAME) $(OPENMP_FLAG_$(F90)) $(LINK_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) \ $(PREFIX) $(LINKERNAME) $(OPENMP_FLAG_$(F90)) $(LINK_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) \
@ -658,8 +715,8 @@ tidy:
@rm -rf *.inst.f90 # for instrumentation @rm -rf *.inst.f90 # for instrumentation
@rm -rf *.pomp.f90 # for instrumentation @rm -rf *.pomp.f90 # for instrumentation
@rm -rf *.pp.f90 # for instrumentation @rm -rf *.pp.f90 # for instrumentation
@rm -rf *.pdb # for instrumnentation @rm -rf *.pdb # for instrumentation
@rm -rf *.opari.inc # for instrumnentation @rm -rf *.opari.inc # for instrumentation
.PHONY: cleanDAMASK .PHONY: cleanDAMASK
cleanDAMASK: cleanDAMASK:

480
code/crystallite.f90
View File

@ -203,11 +203,9 @@ subroutine crystallite_init
tag = '', & tag = '', &
line= '' line= ''
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- crystallite init -+>>>'
write(6,'(/,a)') ' <<<+- crystallite 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
cMax = homogenization_maxNgrains cMax = homogenization_maxNgrains
iMax = mesh_maxNips iMax = mesh_maxNips
@ -520,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, &
@ -584,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
@ -1139,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?
@ -3578,9 +3302,8 @@ logical function crystallite_integrateStress(&
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) &
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress at el ip ipc ',el,ip,ipc write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress at el ip ipc ',el,ip,ipc
endif
#endif #endif
@ -3686,9 +3409,8 @@ logical function crystallite_integrateStress(&
!* calculate plastic velocity gradient and its tangent from constitutive law !* calculate plastic velocity gradient and its tangent from constitutive law
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT3333, dLp_dFi3333, & call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT3333, dLp_dFi3333, &
Tstar_v, Fi_new, ipc, ip, el) Tstar_v, Fi_new, ipc, ip, el)

126
code/damask.core.pyf
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

41
code/debug.f90
View File

@ -104,7 +104,6 @@ contains
subroutine debug_init subroutine debug_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use numerics, only: & use numerics, only: &
worldrank, &
nStress, & nStress, &
nState, & nState, &
nCryst, & nCryst, &
@ -130,47 +129,27 @@ subroutine debug_init
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos
character(len=65536) :: tag, line character(len=65536) :: tag, line
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- debug init -+>>>'
write(6,'(/,a)') ' <<<+- debug 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
if (allocated(debug_StressLoopLpDistribution)) & allocate(debug_StressLoopLpDistribution(nStress+1,2), source=0_pInt)
deallocate(debug_StressLoopLpDistribution) allocate(debug_StressLoopLiDistribution(nStress+1,2), source=0_pInt)
allocate(debug_StressLoopLpDistribution(nStress+1,2)) allocate(debug_StateLoopDistribution(nState+1,2), source=0_pInt)
debug_StressLoopLpDistribution = 0_pInt allocate(debug_CrystalliteLoopDistribution(nCryst+1), source=0_pInt)
if (allocated(debug_StressLoopLiDistribution)) & allocate(debug_MaterialpointStateLoopDistribution(nMPstate), source=0_pInt)
deallocate(debug_StressLoopLiDistribution) allocate(debug_MaterialpointLoopDistribution(nHomog+1), source=0_pInt)
allocate(debug_StressLoopLiDistribution(nStress+1,2))
debug_StressLoopLiDistribution = 0_pInt
if (allocated(debug_StateLoopDistribution)) &
deallocate(debug_StateLoopDistribution)
allocate(debug_StateLoopDistribution(nState+1,2))
debug_StateLoopDistribution = 0_pInt
if (allocated(debug_CrystalliteLoopDistribution)) &
deallocate(debug_CrystalliteLoopDistribution)
allocate(debug_CrystalliteLoopDistribution(nCryst+1))
debug_CrystalliteLoopDistribution = 0_pInt
if (allocated(debug_MaterialpointStateLoopDistribution)) &
deallocate(debug_MaterialpointStateLoopDistribution)
allocate(debug_MaterialpointStateLoopDistribution(nMPstate))
debug_MaterialpointStateLoopDistribution = 0_pInt
if (allocated(debug_MaterialpointLoopDistribution)) &
deallocate(debug_MaterialpointLoopDistribution)
allocate(debug_MaterialpointLoopDistribution(nHomog+1))
debug_MaterialpointLoopDistribution = 0_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! try to open the config file ! try to open the config file
line = '' line = ''
fileExists: if(IO_open_file_stat(FILEUNIT,debug_configFile)) then fileExists: if(IO_open_file_stat(FILEUNIT,debug_configFile)) then
do while (trim(line) /= IO_EOF) ! read thru sections of phase part do while (trim(line) /= IO_EOF) ! read thru sections of phase part
line = IO_read(FILEUNIT) line = IO_read(FILEUNIT)
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag) select case(tag)
case ('element','e','el') case ('element','e','el')
debug_e = IO_intValue(line,chunkPos,2_pInt) debug_e = IO_intValue(line,chunkPos,2_pInt)

18
code/material.f90
View File

@ -499,12 +499,12 @@ subroutine material_init()
allocate(HomogenizationPosition(material_Nhomogenization),source=0_pInt) allocate(HomogenizationPosition(material_Nhomogenization),source=0_pInt)
allocate(CrystallitePosition (material_Nphase), source=0_pInt) allocate(CrystallitePosition (material_Nphase), source=0_pInt)
ElemLoop:do e = 1_pInt,mesh_NcpElems ! loop over elements ElemLoop:do e = 1_pInt,mesh_NcpElems
myHomog = mesh_element(3,e) myHomog = mesh_element(3,e)
IPloop:do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs IPloop:do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
HomogenizationPosition(myHomog) = HomogenizationPosition(myHomog) + 1_pInt HomogenizationPosition(myHomog) = HomogenizationPosition(myHomog) + 1_pInt
mappingHomogenization(1:2,i,e) = [HomogenizationPosition(myHomog),myHomog] mappingHomogenization(1:2,i,e) = [HomogenizationPosition(myHomog),myHomog]
GrainLoop:do g = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) ! loop over grains GrainLoop:do g = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
phase = material_phase(g,i,e) phase = material_phase(g,i,e)
ConstitutivePosition(phase) = ConstitutivePosition(phase)+1_pInt ! not distinguishing between instances of same phase ConstitutivePosition(phase) = ConstitutivePosition(phase)+1_pInt ! not distinguishing between instances of same phase
phaseAt(g,i,e) = phase phaseAt(g,i,e) = phase
@ -1344,10 +1344,10 @@ subroutine material_populateGrains
write(6,'(a32,1x,a32,1x,a6)') 'homogenization_name','microstructure_name','grain#' write(6,'(a32,1x,a32,1x,a6)') 'homogenization_name','microstructure_name','grain#'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
do homog = 1_pInt,material_Nhomogenization ! loop over homogenizations homogenizationLoop: do homog = 1_pInt,material_Nhomogenization
dGrains = homogenization_Ngrains(homog) ! grain number per material point dGrains = homogenization_Ngrains(homog) ! grain number per material point
do micro = 1_pInt,material_Nmicrostructure ! all pairs of homog and micro microstructureLoop: do micro = 1_pInt,material_Nmicrostructure ! all pairs of homog and micro
if (Ngrains(homog,micro) > 0_pInt) then ! an active pair of homog and micro activePair: if (Ngrains(homog,micro) > 0_pInt) then
myNgrains = Ngrains(homog,micro) ! assign short name for total number of grains to populate myNgrains = Ngrains(homog,micro) ! assign short name for total number of grains to populate
myNconstituents = microstructure_Nconstituents(micro) ! assign short name for number of constituents myNconstituents = microstructure_Nconstituents(micro) ! assign short name for number of constituents
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
@ -1578,9 +1578,9 @@ subroutine material_populateGrains
enddo enddo
enddo enddo
endif ! active homog,micro pair endif activePair
enddo enddo microstructureLoop
enddo enddo homogenizationLoop
deallocate(volumeOfGrain) deallocate(volumeOfGrain)
deallocate(phaseOfGrain) deallocate(phaseOfGrain)

58
code/math.f90
View File

@ -158,7 +158,8 @@ module math
math_areaTriangle, & math_areaTriangle, &
math_rotate_forward33, & math_rotate_forward33, &
math_rotate_backward33, & math_rotate_backward33, &
math_rotate_forward3333 math_rotate_forward3333, &
math_limit
private :: & private :: &
math_partition, & math_partition, &
halton, & halton, &
@ -178,7 +179,6 @@ subroutine math_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: tol_math_check use prec, only: tol_math_check
use numerics, only: & use numerics, only: &
worldrank, &
fixedSeed fixedSeed
use IO, only: IO_error, IO_timeStamp use IO, only: IO_error, IO_timeStamp
@ -193,11 +193,9 @@ subroutine math_init
! comment the first random_seed call out, set randSize to 1, and use ifort ! comment the first random_seed call out, set randSize to 1, and use ifort
character(len=64) :: error_msg character(len=64) :: error_msg
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- math init -+>>>'
write(6,'(/,a)') ' <<<+- math 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
call random_seed(size=randSize) call random_seed(size=randSize)
if (allocated(randInit)) deallocate(randInit) if (allocated(randInit)) deallocate(randInit)
@ -216,13 +214,11 @@ subroutine math_init
call random_number(randTest(i)) call random_number(randTest(i))
enddo enddo
mainProcess2: if (worldrank == 0) then write(6,'(a,I2)') ' size of random seed: ', randSize
write(6,*) 'size of random seed: ', randSize do i =1, randSize
do i =1, randSize write(6,'(a,I2,I14)') ' value of random seed: ', i, randInit(i)
write(6,*) 'value of random seed: ', i, randInit(i) enddo
enddo write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest
write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest
endif mainProcess2
call random_seed(put = randInit) call random_seed(put = randInit)
@ -705,7 +701,7 @@ end function math_transpose33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_inv33(A) pure function math_inv33(A)
use prec, only: & use prec, only: &
dNeq dNeq0
implicit none implicit none
real(pReal),dimension(3,3),intent(in) :: A real(pReal),dimension(3,3),intent(in) :: A
@ -718,7 +714,7 @@ pure function math_inv33(A)
DetA = A(1,1) * math_inv33(1,1) + A(1,2) * math_inv33(2,1) + A(1,3) * math_inv33(3,1) DetA = A(1,1) * math_inv33(1,1) + A(1,2) * math_inv33(2,1) + A(1,3) * math_inv33(3,1)
if (dNeq(DetA,0.0_pReal)) then if (dNeq0(DetA)) then
math_inv33(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2) math_inv33(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2)
math_inv33(2,2) = A(1,1) * A(3,3) - A(1,3) * A(3,1) math_inv33(2,2) = A(1,1) * A(3,3) - A(1,3) * A(3,1)
math_inv33(3,2) = -A(1,1) * A(3,2) + A(1,2) * A(3,1) math_inv33(3,2) = -A(1,1) * A(3,2) + A(1,2) * A(3,1)
@ -743,7 +739,7 @@ end function math_inv33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine math_invert33(A, InvA, DetA, error) pure subroutine math_invert33(A, InvA, DetA, error)
use prec, only: & use prec, only: &
dEq dEq0
implicit none implicit none
logical, intent(out) :: error logical, intent(out) :: error
@ -757,7 +753,7 @@ pure subroutine math_invert33(A, InvA, DetA, error)
DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1) DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1)
if (dEq(DetA,0.0_pReal)) then if (dEq0(DetA)) then
InvA = 0.0_pReal InvA = 0.0_pReal
error = .true. error = .true.
else else
@ -1081,7 +1077,7 @@ pure function math_Plain99to3333(m99)
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain99to3333(mapPlain(1,i),mapPlain(2,i),& forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain99to3333(mapPlain(1,i),mapPlain(2,i),&
mapPlain(1,j),mapPlain(2,j)) = m99(i,j) mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
end function math_Plain99to3333 end function math_Plain99to3333
@ -1193,10 +1189,10 @@ function math_qRand()
real(pReal), dimension(3) :: rnd real(pReal), dimension(3) :: rnd
call halton(3_pInt,rnd) call halton(3_pInt,rnd)
math_qRand(1) = cos(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)) math_qRand = [cos(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)), &
math_qRand(2) = sin(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)) sin(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
math_qRand(3) = cos(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)) cos(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
math_qRand(4) = sin(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)) sin(2.0_pReal*PI*rnd(1))*sqrt(rnd(3))]
end function math_qRand end function math_qRand
@ -1263,7 +1259,7 @@ end function math_qNorm
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_qInv(Q) pure function math_qInv(Q)
use prec, only: & use prec, only: &
dNeq dNeq0
implicit none implicit none
real(pReal), dimension(4), intent(in) :: Q real(pReal), dimension(4), intent(in) :: Q
@ -1273,7 +1269,7 @@ pure function math_qInv(Q)
math_qInv = 0.0_pReal math_qInv = 0.0_pReal
squareNorm = math_qDot(Q,Q) squareNorm = math_qDot(Q,Q)
if (dNeq(squareNorm,0.0_pReal)) math_qInv = math_qConj(Q) / squareNorm if (dNeq0(squareNorm)) math_qInv = math_qConj(Q) / squareNorm
end function math_qInv end function math_qInv
@ -1474,7 +1470,7 @@ pure function math_axisAngleToR(axis,omega)
real(pReal), dimension(3) :: axisNrm real(pReal), dimension(3) :: axisNrm
real(pReal) :: norm,s,c,c1 real(pReal) :: norm,s,c,c1
norm = sqrt(math_mul3x3(axis,axis)) norm = norm2(axis)
if (norm > 1.0e-8_pReal) then ! non-zero rotation if (norm > 1.0e-8_pReal) then ! non-zero rotation
axisNrm = axis/norm ! normalize axis to be sure axisNrm = axis/norm ! normalize axis to be sure
@ -1580,7 +1576,7 @@ pure function math_qToR(q)
S = reshape( [0.0_pReal, -q(4), q(3), & S = reshape( [0.0_pReal, -q(4), q(3), &
q(4), 0.0_pReal, -q(2), & q(4), 0.0_pReal, -q(2), &
-q(3), q(2), 0.0_pReal],[3,3]) ! notation is transposed -q(3), q(2), 0.0_pReal],[3,3]) ! notation is transposed
math_qToR = (2.0_pReal * q(1)*q(1) - 1.0_pReal) * math_I3 & math_qToR = (2.0_pReal * q(1)*q(1) - 1.0_pReal) * math_I3 &
+ 2.0_pReal * T - 2.0_pReal * q(1) * S + 2.0_pReal * T - 2.0_pReal * q(1) * S
@ -1604,17 +1600,17 @@ pure function math_qToEuler(qPassive)
q = math_qConj(qPassive) ! convert to active rotation, since formulas are defined for active rotations q = math_qConj(qPassive) ! convert to active rotation, since formulas are defined for active rotations
math_qToEuler(2) = acos(1.0_pReal-2.0_pReal*(q(2)*q(2)+q(3)*q(3))) math_qToEuler(2) = acos(1.0_pReal-2.0_pReal*(q(2)**2+q(3)**2))
if (abs(math_qToEuler(2)) < 1.0e-6_pReal) then if (abs(math_qToEuler(2)) < 1.0e-6_pReal) then
math_qToEuler(1) = sign(2.0_pReal*acos(math_limit(q(1),-1.0_pReal, 1.0_pReal)),q(4)) math_qToEuler(1) = sign(2.0_pReal*acos(math_limit(q(1),-1.0_pReal, 1.0_pReal)),q(4))
math_qToEuler(3) = 0.0_pReal math_qToEuler(3) = 0.0_pReal
else else
math_qToEuler(1) = atan2(q(1)*q(3)+q(2)*q(4), q(1)*q(2)-q(3)*q(4)) math_qToEuler(1) = atan2(+q(1)*q(3)+q(2)*q(4), q(1)*q(2)-q(3)*q(4))
math_qToEuler(3) = atan2(-q(1)*q(3)+q(2)*q(4), q(1)*q(2)+q(3)*q(4)) math_qToEuler(3) = atan2(-q(1)*q(3)+q(2)*q(4), q(1)*q(2)+q(3)*q(4))
endif endif
math_qToEuler = merge(math_qToEuler + [2.0_pReal*PI, PI, 2.0_pReal*PI], & ! ensure correct range math_qToEuler = merge(math_qToEuler + [2.0_pReal*PI, PI, 2.0_pReal*PI], & ! ensure correct range
math_qToEuler, math_qToEuler<0.0_pReal) math_qToEuler, math_qToEuler<0.0_pReal)
end function math_qToEuler end function math_qToEuler
@ -2078,7 +2074,7 @@ end function math_eigenvectorBasisSym33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function math_rotationalPart33(m) function math_rotationalPart33(m)
use prec, only: & use prec, only: &
dEq dEq0
use IO, only: & use IO, only: &
IO_warning IO_warning
@ -2090,7 +2086,7 @@ function math_rotationalPart33(m)
U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m)) U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m))
Uinv = math_inv33(U) Uinv = math_inv33(U)
inversionFailed: if (all(dEq(Uinv,0.0_pReal))) then inversionFailed: if (all(dEq0(Uinv))) then
math_rotationalPart33 = math_I3 math_rotationalPart33 = math_I3
call IO_warning(650_pInt) call IO_warning(650_pInt)
else inversionFailed else inversionFailed

13
code/mesh.f90
View File

@ -481,6 +481,7 @@ subroutine mesh_init(ip,el)
#endif #endif
#ifdef Spectral #ifdef Spectral
IO_open_file, & IO_open_file, &
IO_error, &
#else #else
IO_open_InputFile, & IO_open_InputFile, &
#endif #endif
@ -507,7 +508,8 @@ subroutine mesh_init(ip,el)
implicit none implicit none
#ifdef Spectral #ifdef Spectral
integer(C_INTPTR_T) :: gridMPI(3), alloc_local, local_K, local_K_offset integer(C_INTPTR_T) :: devNull, local_K, local_K_offset
integer :: ierr, worldsize
#endif #endif
integer(pInt), parameter :: FILEUNIT = 222_pInt integer(pInt), parameter :: FILEUNIT = 222_pInt
integer(pInt), intent(in) :: el, ip integer(pInt), intent(in) :: el, ip
@ -547,10 +549,13 @@ subroutine mesh_init(ip,el)
call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file... call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file...
if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6) if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6)
grid = mesh_spectral_getGrid(fileUnit) grid = mesh_spectral_getGrid(fileUnit)
call MPI_comm_size(MPI_COMM_WORLD, worldsize, ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_comm_size')
if(worldsize>grid(3)) call IO_error(894_pInt, ext_msg='number of processes exceeds grid(3)')
geomSize = mesh_spectral_getSize(fileUnit) geomSize = mesh_spectral_getSize(fileUnit)
gridMPI = int(grid,C_INTPTR_T) devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T),int(grid(2),C_INTPTR_T),&
alloc_local = fftw_mpi_local_size_3d(gridMPI(3), gridMPI(2), gridMPI(1)/2 +1, & int(grid(1),C_INTPTR_T)/2+1,MPI_COMM_WORLD,local_K,local_K_offset)
MPI_COMM_WORLD, local_K, local_K_offset)
grid3 = int(local_K,pInt) grid3 = int(local_K,pInt)
grid3Offset = int(local_K_offset,pInt) grid3Offset = int(local_K_offset,pInt)
size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal)

4
code/numerics.f90
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

12
code/plastic_disloUCLA.f90
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)

12
code/plastic_dislotwin.f90
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)

12
code/plastic_isotropic.f90
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)

9
code/plastic_none.f90
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))

14
code/plastic_nonlocal.f90
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
@ -1119,7 +1115,7 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances),
initializeInstances: do phase = 1_pInt, size(phase_plasticity) initializeInstances: do phase = 1_pInt, size(phase_plasticity)
NofMyPhase=count(material_phase==phase) NofMyPhase=count(material_phase==phase)
myPhase2: if (phase_plasticity(phase) == PLASTICITY_NONLOCAL_ID .and. NofMyPhase/=0) then myPhase2: if (phase_plasticity(phase) == PLASTICITY_NONLOCAL_ID) then
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
!*** Inverse lookup of my slip system family and the slip system in lattice !*** Inverse lookup of my slip system family and the slip system in lattice
@ -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)

12
code/plastic_phenoplus.f90
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)

18
code/plastic_phenopowerlaw.f90
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)
@ -1160,9 +1152,9 @@ function plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
enddo enddo
plastic_phenopowerlaw_postResults(c+j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* & plastic_phenopowerlaw_postResults(c+j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos)/plasticState(ph)%state(j,of))**plastic_phenopowerlaw_n_slip(instance) & ((abs(tau_slip_pos)/plasticState(ph)%state(j,of))**plastic_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/plasticState(ph)%state(j,of))**plastic_phenopowerlaw_n_slip(instance))& *sign(1.0_pReal,tau_slip_pos) &
*sign(1.0_pReal,tau_slip_pos) +(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) &
*sign(1.0_pReal,tau_slip_neg))
enddo slipSystems1 enddo slipSystems1
enddo slipFamilies1 enddo slipFamilies1
c = c + nSlip c = c + nSlip

12
code/plastic_titanmod.f90
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)

39
code/prec.f90
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
@ -83,7 +81,7 @@ module prec
nTwin = 0_pInt, & nTwin = 0_pInt, &
nTrans = 0_pInt nTrans = 0_pInt
logical :: & logical :: &
nonlocal = .false. !< absolute tolerance for state integration nonlocal = .false.
real(pReal), pointer, dimension(:,:), contiguous :: & real(pReal), pointer, dimension(:,:), contiguous :: &
slipRate, & !< slip rate slipRate, & !< slip rate
accumulatedSlip !< accumulated plastic slip accumulatedSlip !< accumulated plastic slip
@ -115,8 +113,10 @@ module prec
prec_init, & prec_init, &
prec_isNaN, & prec_isNaN, &
dEq, & dEq, &
dEq0, &
cEq, & cEq, &
dNeq, & dNeq, &
dNeq0, &
cNeq cNeq
contains contains
@ -199,6 +199,39 @@ logical elemental pure function dNeq(a,b,tol)
dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
end function dNeq end function dNeq
!--------------------------------------------------------------------------------------------------
!> @brief equality to 0comparison for float with double precision
! replaces "==0" but for certain (relative) tolerance. Counterpart to dNeq0
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
!--------------------------------------------------------------------------------------------------
logical elemental pure function dEq0(a,tol)
implicit none
real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
dEq0 = merge(.True., .False.,abs(a) <= merge(tol,eps,present(tol))*abs(a))
end function dEq0
!--------------------------------------------------------------------------------------------------
!> @brief inequality comparison to 0 for float with double precision
! replaces "!=0" but for certain (relative) tolerance. Counterpart to dEq0
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
!--------------------------------------------------------------------------------------------------
logical elemental pure function dNeq0(a,tol)
implicit none
real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: tol
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
dNeq0 = merge(.False., .True.,abs(a) <= merge(tol,eps,present(tol))*abs(a))
end function dNeq0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief equality comparison for complex with double precision !> @brief equality comparison for complex with double precision
! replaces "==" but for certain (relative) tolerance. Counterpart to cNeq ! replaces "==" but for certain (relative) tolerance. Counterpart to cNeq

12
code/source_damage_anisoBrittle.f90
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)

12
code/source_damage_anisoDuctile.f90
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)

12
code/source_damage_isoBrittle.f90
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)

12
code/source_damage_isoDuctile.f90
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)

12
code/source_thermal_dissipation.f90
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)

12
code/source_thermal_externalheat.f90
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)

12
code/source_vacancy_irradiation.f90
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)

12
code/source_vacancy_phenoplasticity.f90
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)

12
code/source_vacancy_thermalfluc.f90
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)

11
code/spectral_interface.f90
View File

@ -57,7 +57,8 @@ subroutine DAMASK_interface_init()
integer :: & integer :: &
i, & i, &
threadLevel, & threadLevel, &
worldrank = 0 worldrank = 0, &
worldsize = 0
integer, allocatable, dimension(:) :: & integer, allocatable, dimension(:) :: &
chunkPos chunkPos
integer, dimension(8) :: & integer, dimension(8) :: &
@ -66,6 +67,7 @@ subroutine DAMASK_interface_init()
external :: & external :: &
quit,& quit,&
MPI_Comm_rank,& MPI_Comm_rank,&
MPI_Comm_size,&
PETScInitialize, & PETScInitialize, &
MPI_Init_Thread, & MPI_Init_Thread, &
MPI_abort MPI_abort
@ -77,17 +79,17 @@ subroutine DAMASK_interface_init()
#ifdef _OPENMP #ifdef _OPENMP
call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr) ! in case of OpenMP, don't rely on PETScInitialize doing MPI init call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr) ! in case of OpenMP, don't rely on PETScInitialize doing MPI init
if (threadLevel<MPI_THREAD_FUNNELED) then if (threadLevel<MPI_THREAD_FUNNELED) then
write(6,'(a)') 'MPI library does not support OpenMP' write(6,'(a)') ' MPI library does not support OpenMP'
call quit(1_pInt) call quit(1_pInt)
endif endif
#endif #endif
call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code
CHKERRQ(ierr) ! this is a macro definition, it is case sensitive CHKERRQ(ierr) ! this is a macro definition, it is case sensitive
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr) call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
call MPI_Comm_size(MPI_COMM_WORLD, worldsize, ierr);CHKERRQ(ierr)
mainProcess: if (worldrank == 0) then mainProcess: if (worldrank == 0) then
if (output_unit /= 6) then if (output_unit /= 6) then
write(output_unit,'(a)') 'STDOUT != 6' write(output_unit,'(a)') ' STDOUT != 6'
call quit(1_pInt) call quit(1_pInt)
endif endif
else mainProcess else mainProcess
@ -104,6 +106,7 @@ subroutine DAMASK_interface_init()
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',& write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',&
dateAndTime(6),':',& dateAndTime(6),':',&
dateAndTime(7) dateAndTime(7)
write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
#include "compilation_info.f90" #include "compilation_info.f90"

2
code/spectral_thermal.f90
View File

@ -245,7 +245,7 @@ type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_ol
if (worldrank == 0) then if (worldrank == 0) then
if (spectral_thermal_solution%converged) & if (spectral_thermal_solution%converged) &
write(6,'(/,a)') ' ... thermal conduction converged ..................................' write(6,'(/,a)') ' ... thermal conduction converged ..................................'
write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature = ',& write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',&
minTemperature, maxTemperature, stagNorm minTemperature, maxTemperature, stagNorm
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)

104
code/spectral_utilities.f90
View File

@ -172,8 +172,7 @@ subroutine utilities_init()
memory_efficient, & memory_efficient, &
petsc_defaultOptions, & petsc_defaultOptions, &
petsc_options, & petsc_options, &
divergence_correction, & divergence_correction
worldrank
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_SPECTRAL, & debug_SPECTRAL, &
@ -212,11 +211,9 @@ subroutine utilities_init()
vecSize = 3_C_INTPTR_T, & vecSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T tensorSize = 9_C_INTPTR_T
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>'
write(6,'(/,a)') ' <<<+- spectral_utilities 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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set debugging parameters ! set debugging parameters
@ -224,11 +221,11 @@ subroutine utilities_init()
debugRotation = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION) /= 0 debugRotation = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION) /= 0
debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0 debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0
if(debugPETSc .and. worldrank == 0_pInt) write(6,'(3(/,a),/)') & if(debugPETSc) write(6,'(3(/,a),/)') &
' Initializing PETSc with debug options: ', & ' Initializing PETSc with debug options: ', &
trim(PETScDebug), & trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.config ' ' add more using the PETSc_Options keyword in numerics.config '; flush(6)
flush(6)
call PetscOptionsClear(ierr); CHKERRQ(ierr) call PetscOptionsClear(ierr); CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(trim(PETSCDEBUG),ierr); CHKERRQ(ierr) if(debugPETSc) call PetscOptionsInsertString(trim(PETSCDEBUG),ierr); CHKERRQ(ierr)
call PetscOptionsInsertString(trim(petsc_defaultOptions),ierr); CHKERRQ(ierr) call PetscOptionsInsertString(trim(petsc_defaultOptions),ierr); CHKERRQ(ierr)
@ -237,10 +234,8 @@ subroutine utilities_init()
grid1Red = grid(1)/2_pInt + 1_pInt grid1Red = grid(1)/2_pInt + 1_pInt
wgt = 1.0/real(product(grid),pReal) wgt = 1.0/real(product(grid),pReal)
if (worldrank == 0) then write(6,'(a,3(i12 ))') ' grid a b c: ', grid
write(6,'(a,3(i12 ))') ' grid a b c: ', grid write(6,'(a,3(es12.5))') ' size x y z: ', geomSize
write(6,'(a,3(es12.5))') ' size x y z: ', geomSize
endif
select case (spectral_derivative) select case (spectral_derivative)
case ('continuous') ! default, no weighting case ('continuous') ! default, no weighting
@ -342,8 +337,7 @@ subroutine utilities_init()
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(0_pInt,ext_msg='Fortran to C') ! check for correct precision in C if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(0_pInt,ext_msg='Fortran to C') ! check for correct precision in C
call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation
if (debugGeneral .and. worldrank == 0_pInt) write(6,'(/,a)') ' FFTW initialized' if (debugGeneral) write(6,'(/,a)') ' FFTW initialized'; flush(6)
flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
@ -527,8 +521,6 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
use math, only: & use math, only: &
math_det33, & math_det33, &
math_invert math_invert
use numerics, only: &
worldrank
use mesh, only: & use mesh, only: &
grid3, & grid3, &
grid, & grid, &
@ -545,10 +537,8 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
logical :: err logical :: err
if (worldrank == 0_pInt) then write(6,'(/,a)') ' ... doing gamma convolution ...............................................'
write(6,'(/,a)') ' ... doing gamma convolution ...............................................' flush(6)
flush(6)
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium) ! do the actual spectral method calculation (mechanical equilibrium)
@ -624,8 +614,6 @@ end subroutine utilities_fourierGreenConvolution
real(pReal) function utilities_divergenceRMS() real(pReal) function utilities_divergenceRMS()
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: &
worldrank
use mesh, only: & use mesh, only: &
geomSize, & geomSize, &
grid, & grid, &
@ -638,10 +626,9 @@ real(pReal) function utilities_divergenceRMS()
external :: & external :: &
MPI_Allreduce MPI_Allreduce
if (worldrank == 0_pInt) then write(6,'(/,a)') ' ... calculating divergence ................................................'
write(6,'(/,a)') ' ... calculating divergence ................................................' flush(6)
flush(6)
endif
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -680,8 +667,6 @@ end function utilities_divergenceRMS
real(pReal) function utilities_curlRMS() real(pReal) function utilities_curlRMS()
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: &
worldrank
use mesh, only: & use mesh, only: &
geomSize, & geomSize, &
grid, & grid, &
@ -693,13 +678,11 @@ real(pReal) function utilities_curlRMS()
complex(pReal), dimension(3) :: rescaledGeom complex(pReal), dimension(3) :: rescaledGeom
external :: & external :: &
MPI_Reduce, &
MPI_Allreduce MPI_Allreduce
if (worldrank == 0_pInt) then write(6,'(/,a)') ' ... calculating curl ......................................................'
write(6,'(/,a)') ' ... calculating curl ......................................................' flush(6)
flush(6)
endif
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -757,8 +740,6 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
prec_isNaN prec_isNaN
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: &
worldrank
use math, only: & use math, only: &
math_Plain3333to99, & math_Plain3333to99, &
math_plain99to3333, & math_plain99to3333, &
@ -790,7 +771,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal) allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal)
temp99_Real = math_Plain3333to99(math_rotate_forward3333(C,rot_BC)) temp99_Real = math_Plain3333to99(math_rotate_forward3333(C,rot_BC))
if(debugGeneral .and. worldrank == 0_pInt) then if(debugGeneral) then
write(6,'(/,a)') ' ... updating masked compliance ............................................' write(6,'(/,a)') ' ... updating masked compliance ............................................'
write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',&
transpose(temp99_Real)/1.e9_pReal transpose(temp99_Real)/1.e9_pReal
@ -831,7 +812,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
if(m/=n .and. abs(sTimesC(m,n)) > (0.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! off diagonal elements of S*C should be 0 if(m/=n .and. abs(sTimesC(m,n)) > (0.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! off diagonal elements of S*C should be 0
enddo enddo
enddo enddo
if((debugGeneral .or. errmatinv) .and. (worldrank == 0_pInt)) then ! report if(debugGeneral .or. errmatinv) then
write(formatString, '(I16.16)') size_reduced write(formatString, '(I16.16)') size_reduced
formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
write(6,trim(formatString),advance='no') ' C * S (load) ', & write(6,trim(formatString),advance='no') ' C * S (load) ', &
@ -845,7 +826,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
else else
temp99_real = 0.0_pReal temp99_real = 0.0_pReal
endif endif
if(debugGeneral .and. worldrank == 0_pInt) & ! report if(debugGeneral) &
write(6,'(/,a,/,9(9(2x,f12.7,1x)/),/)',advance='no') ' Masked Compliance (load) * GPa =', & write(6,'(/,a,/,9(9(2x,f12.7,1x)/),/)',advance='no') ' Masked Compliance (load) * GPa =', &
transpose(temp99_Real*1.e9_pReal) transpose(temp99_Real*1.e9_pReal)
flush(6) flush(6)
@ -938,15 +919,11 @@ end subroutine utilities_fourierTensorDivergence
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC) P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC)
use prec, only: &
dNeq
use IO, only: & use IO, only: &
IO_error IO_error
use debug, only: & use debug, only: &
debug_reset, & debug_reset, &
debug_info debug_info
use numerics, only: &
worldrank
use math, only: & use math, only: &
math_transpose33, & math_transpose33, &
math_rotate_forward33, & math_rotate_forward33, &
@ -974,7 +951,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress
logical :: & logical :: &
age age
@ -984,26 +961,18 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
external :: & write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
MPI_Reduce, & flush(6)
MPI_Allreduce
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
flush(6)
endif
age = .False. age = .False.
if (forwardData) then ! aging results if (forwardData) then ! aging results
age = .True. age = .True.
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3])
endif endif
if (cutBack) then ! restore saved variables if (cutBack) age = .False. ! restore saved variables
age = .False.
endif
materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3])
call debug_reset() call debug_reset() ! this has no effect on rank >0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate bounds of det(F) and report ! calculate bounds of det(F) and report
@ -1015,15 +984,10 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
defgradDetMax = max(defgradDetMax,defgradDet) defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet) defgradDetMin = min(defgradDetMin,defgradDet)
end do end do
call MPI_reduce(MPI_IN_PLACE,defgradDetMax,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max') write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
call MPI_reduce(MPI_IN_PLACE,defgradDetMin,1,MPI_DOUBLE,MPI_MIN,0,PETSC_COMM_WORLD,ierr) write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min') flush(6)
if (worldrank == 0_pInt) then
write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
flush(6)
endif
endif endif
call CPFEM_general(age,timeinc) call CPFEM_general(age,timeinc)
@ -1053,7 +1017,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
call debug_info() call debug_info() ! this has no effect on rank >0
restartWrite = .false. ! reset restartWrite status restartWrite = .false. ! reset restartWrite status
cutBack = .false. ! reset cutBack status cutBack = .false. ! reset cutBack status
@ -1061,15 +1025,13 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (debugRotation .and. worldrank == 0_pInt) & if (debugRotation) &
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
math_transpose33(P_av)*1.e-6_pReal math_transpose33(P_av)*1.e-6_pReal
P_av = math_rotate_forward33(P_av,rotation_BC) P_av = math_rotate_forward33(P_av,rotation_BC)
if (worldrank == 0_pInt) then write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
math_transpose33(P_av)*1.e-6_pReal math_transpose33(P_av)*1.e-6_pReal
flush(6) flush(6)
endif
end subroutine utilities_constitutiveResponse end subroutine utilities_constitutiveResponse

35
config/Phase_Phenopowerlaw_BCC-Ferrite.config
View File

@ -1,35 +0,0 @@
# Tasan et.al. 2015 Acta Materalia
# Tasan et.al. 2015 International Journal of Plasticity
# Diehl et.al. 2015 Meccanica
[BCC_Ferrite]
elasticity hooke
plasticity phenopowerlaw
lattice_structure bcc
Nslip 12 12 # per family
Ntwin 0 # per family
c11 233.3e9
c12 135.5e9
c44 118.0e9
gdot0_slip 0.001
n_slip 20
tau0_slip 95.e6 97.e6 0 0 # per family, optimization long simplex 109
tausat_slip 222.e6 412.7e6 0 0 # per family, optimization long simplex 109
gdot0_twin 0.001
n_twin 20
tau0_twin 31.0e6 # per family
s_pr 0 # push-up factor for slip saturation due to twinning
twin_b 0
twin_c 0
twin_d 0
twin_e 0
h0_slipslip 1000.0e6 # opti
h0_twinslip 0
h0_twintwin 0
interaction_slipslip 1 1 1.4 1.4 1.4 1.4
interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
w0_slip 2.0 # opti
atol_resistance 1
(output) totalshear

138
configure vendored
View File

@ -1,138 +0,0 @@
#!/usr/bin/env python
import os,re,sys,string,subprocess,shutil
from optparse import OptionParser, Option
# -----------------------------
class extendableOption(Option):
# -----------------------------
# used for definition of new option parser action 'extend', which enables to take multiple option arguments
# taken from online tutorial http://docs.python.org/library/optparse.html
ACTIONS = Option.ACTIONS + ("extend",)
STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
def take_action(self, action, dest, opt, value, values, parser):
if action == "extend":
lvalue = value.split(",")
values.ensure_value(dest, []).extend(lvalue)
else:
Option.take_action(self, action, dest, opt, value, values, parser)
########################################################
# MAIN
########################################################
parser = OptionParser(option_class=extendableOption, usage='%prog options', description = """
Configures the compilation and installation of DAMASK
""")
defaults={'DAMASK_BIN':'depending on access rights',
'FFTW_ROOT':'/usr',
'MSC_ROOT' :'/msc',
'DAMASK_NUM_THREADS':4,
'MARC_VERSION':'2015',
'spectralOptions':{},
}
#--- if local config file exists, read, otherwise assume global config file ------------------------
configFile = os.path.join(os.getenv('HOME'),'.damask/damask.conf') \
if os.path.isfile(os.path.join(os.getenv('HOME'),'.damask/damask.conf')) \
else '/etc/damask.conf'
#--- set default values according to read in values ------------------------------------------------
try:
with open(configFile,'r') as f:
print('\n<<<<< reading default values from %s\n'%configFile)
for line in f:
line = line.strip()
if line.startswith('#') or line == '':
pass
[key,value] = (re.split('[= ]',line)+[None,None])[:2]
if key == 'DAMASK_NUM_THREADS':
defaults['DAMASK_NUM_THREADS'] = int(value)
if key == 'DAMASK_BIN':
defaults['DAMASK_BIN'] = value
except IOError:
pass
parser.add_option('--prefix', dest='prefix', metavar='string',
help='location of (links to) DAMASK executables [%default]')
parser.add_option('--with-OMP-threads','--with-omp-threads',
dest='threads', type='int', metavar='int',
help='number of openMP threads [%default]')
parser.add_option('--with-spectral-options', dest='spectraloptions', action='extend', metavar='<string LIST>',
help='options for compilation of spectral solver')
parser.set_defaults(prefix = defaults['DAMASK_BIN'])
parser.set_defaults(mscRoot = defaults['MSC_ROOT'])
parser.set_defaults(marcVersion = defaults['MARC_VERSION'])
parser.set_defaults(threads = defaults['DAMASK_NUM_THREADS'])
parser.set_defaults(spectraloptions = [])
(options,filenames) = parser.parse_args()
#--- read config file if present to keep comments and order ---------------------------------------
output = []
try:
with open(configFile,'r') as f:
for line in f:
line = line.strip()
items = re.split('[= ]',line)
if (not line or items[0].startswith('#')):
pass
if items[0] == 'DAMASK_BIN':
line = '%s=%s'%(items[0],options.prefix)
options.prefix ='depending on access rights'
if items[0] == 'MSC_ROOT':
line = '%s=%s'%(items[0],options.mscRoot)
options.mscRoot =''
if items[0] == 'MARC_VERSION':
line = '%s=%s'%(items[0],options.marcVersion)
options.marcVersion =''
if items[0] == 'DAMASK_NUM_THREADS':
line = '%s=%s'%(items[0],options.threads)
options.threads =''
for spectralOption in options.spectraloptions:
[key,value] = re.split('[= ]',spectralOption)[0:2]
if key == items[0]:
line = '%s=%s'%(items[0],value)
options.spectraloptions.remove(spectralOption)
output.append(line)
except IOError:
pass
#--- write remaining options --------------------------------------------------------------------------
for opt, value in options.__dict__.items():
if opt == 'prefix' and value != 'depending on access rights':
output.append('DAMASK_BIN=%s'%value)
if opt == 'mscRoot' and value != '':
output.append('MSC_ROOT=%s'%value)
if opt == 'marcVersion' and value != '':
output.append('MARC_VERSION=%s'%value)
if opt == 'threads' and value != '':
output.append('DAMASK_NUM_THREADS=%s'%value)
for spectralOption in options.spectraloptions:
output.append(spectralOption)
#--- decide where do save the data -------------------------------------------------------------------
configDir = '/etc' if os.access('/etc/', os.W_OK) \
else os.path.join(os.getenv('HOME'),'.damask') # use system-wide config if possible
configFileNew = os.path.join(configDir,'damask.conf')
if not os.path.isdir(configDir):
os.mkdir(configDir)
print('\n>>>>> writing values to %s\n'%configFileNew)
with open(configFileNew,'w') as f:
for line in output:
print(line)
f.write(line+'\n')

1
examples/AbaqusStandard/numerics.config
View File

@ -1,2 +1 @@
fixed_seed 1697667030 fixed_seed 1697667030
analyticJaco 1

0
config/Crystallite_All.config → examples/ConfigFiles/Crystallite_All.config
View File

0
config/Crystallite_None.config → examples/ConfigFiles/Crystallite_None.config
View File

0
config/Crystallite_aLittleSomething.config → examples/ConfigFiles/Crystallite_aLittleSomething.config
View File

0
config/Homogenization_Damage_NonLocal.config → examples/ConfigFiles/Homogenization_Damage_NonLocal.config
View File

0
config/Homogenization_HydrogenFlux_CahnHilliard.config → examples/ConfigFiles/Homogenization_HydrogenFlux_CahnHilliard.config
View File

0
config/Homogenization_Isostrain_Parallel3.config → examples/ConfigFiles/Homogenization_Isostrain_Parallel3.config
View File

0
config/Homogenization_Isostrain_SX.config → examples/ConfigFiles/Homogenization_Isostrain_SX.config
View File

0
config/Homogenization_Isostrain_Taylor2.config → examples/ConfigFiles/Homogenization_Isostrain_Taylor2.config
View File

0
config/Homogenization_None_Dummy.config → examples/ConfigFiles/Homogenization_None_Dummy.config
View File

0
config/Homogenization_Porosity_PhaseField.config → examples/ConfigFiles/Homogenization_Porosity_PhaseField.config
View File

0
config/Homogenization_RGC_8Grains.config → examples/ConfigFiles/Homogenization_RGC_8Grains.config
View File

0
config/Homogenization_Thermal_Conduction.config → examples/ConfigFiles/Homogenization_Thermal_Conduction.config
View File

0
config/Homogenization_VacancyFlux_CahnHilliard.config → examples/ConfigFiles/Homogenization_VacancyFlux_CahnHilliard.config
View File

0
config/Homogenization_multiField.config → examples/ConfigFiles/Homogenization_multiField.config
View File

0
config/Kinematics_Hydrogen_Strain.config → examples/ConfigFiles/Kinematics_Hydrogen_Strain.config
View File

0
config/Kinematics_Thermal_Expansion.config → examples/ConfigFiles/Kinematics_Thermal_Expansion.config
View File

0
config/Kinematics_Vacancy_Strain.config → examples/ConfigFiles/Kinematics_Vacancy_Strain.config
View File

0
config/Microstructure_DP_Steel.config → examples/ConfigFiles/Microstructure_DP_Steel.config
View File

0
config/Microstructure_ElementHomogeneous.config → examples/ConfigFiles/Microstructure_ElementHomogeneous.config
View File

0
config/Phase_Damage.config → examples/ConfigFiles/Phase_Damage.config
View File

0
config/Phase_DisloUCLA_Tungsten.config → examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config
View File

0
config/Phase_Dislotwin_TWIP-Steel-FeMnC.config → examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config
View File

0
config/Phase_Dislotwin_Tungsten.config → examples/ConfigFiles/Phase_Dislotwin_Tungsten.config
View File

0
config/Phase_Hydrogen.config → examples/ConfigFiles/Phase_Hydrogen.config
View File

0
config/Phase_Isotropic_AluminumIsotropic.config → examples/ConfigFiles/Phase_Isotropic_AluminumIsotropic.config
View File

0
config/Phase_None_IsotropicVolumePreservation.config → examples/ConfigFiles/Phase_None_IsotropicVolumePreservation.config
View File

0
config/Phase_None_Orthorombic.config → examples/ConfigFiles/Phase_None_Orthorombic.config
View File

0
config/Phase_Nonlocal_Aluminum.config → examples/ConfigFiles/Phase_Nonlocal_Aluminum.config
View File

94
examples/ConfigFiles/Phase_Nonlocal_Nickel.config Executable file
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

0
config/Phase_Phenopowerlaw_Aluminum.config → examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config
View File

25
examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config Normal file
View File

@ -0,0 +1,25 @@
# Tasan et.al. 2015 Acta Materalia
# Tasan et.al. 2015 International Journal of Plasticity
# Diehl et.al. 2015 Meccanica
[BCC-Ferrite]
elasticity hooke
plasticity phenopowerlaw
lattice_structure bcc
Nslip 12 12 # per family
Ntwin 0 # per family
c11 233.3e9
c12 135.5e9
c44 118.0e9
gdot0_slip 0.001
n_slip 20
tau0_slip 95.e6 97.e6 # per family, optimization long simplex 109
tausat_slip 222.e6 412.7e6 # per family, optimization long simplex 109
h0_slipslip 1000.0e6
interaction_slipslip 1 1 1.4 1.4 1.4 1.4
interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
w0_slip 2.0
(output) totalshear

25
config/Phase_Phenopowerlaw_BCC-Martensite.config → examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Martensite.config
View File

@ -1,34 +1,25 @@
# Tasan et.al. 2015 Acta Materalia # Tasan et.al. 2015 Acta Materalia
# Tasan et.al. 2015 International Journal of Plasticity # Tasan et.al. 2015 International Journal of Plasticity
# Diehl et.al. 2015 Meccanica # Diehl et.al. 2015 Meccanica
[BCC_Martensite] [BCC-Martensite]
plasticity phenopowerlaw
elasticity hooke elasticity hooke
plasticity phenopowerlaw
lattice_structure bcc lattice_structure bcc
Nslip 12 12 # per family Nslip 12 12 # per family
Ntwin 0 # per family Ntwin 0 # per family
c11 417.4e9 c11 417.4e9
c12 242.4e9 c12 242.4e9
c44 211.1e9 c44 211.1e9
gdot0_slip 0.001 gdot0_slip 0.001
n_slip 20 n_slip 20
tau0_slip 405.8e6 456.7e6 0 0 # per family tau0_slip 405.8e6 456.7e6 # per family
tausat_slip 872.9e6 971.2e6 0 0 # per family tausat_slip 872.9e6 971.2e6 # per family
gdot0_twin 0.001
n_twin 20
tau0_twin 31.0e6 # per family
s_pr 0 # push-up factor for slip saturation due to twinning
twin_b 0
twin_c 0
twin_d 0
twin_e 0
h0_slipslip 563.0e9 h0_slipslip 563.0e9
h0_twinslip 0
h0_twintwin 0
interaction_slipslip 1 1 1.4 1.4 1.4 1.4 interaction_slipslip 1 1 1.4 1.4 1.4 1.4
interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
w0_slip 2.0 w0_slip 2.0
atol_resistance 1
(output) totalshear (output) totalshear

0
config/Phase_Phenopowerlaw_Gold.config → examples/ConfigFiles/Phase_Phenopowerlaw_Gold.config
View File

0
config/Phase_Phenopowerlaw_Magnesium.config → examples/ConfigFiles/Phase_Phenopowerlaw_Magnesium.config
View File

0
config/Phase_Phenopowerlaw_cpTi-alpha.config → examples/ConfigFiles/Phase_Phenopowerlaw_cpTi-alpha.config
View File

0
config/Phase_Phenopowerlaw_multiField.config → examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config
View File

0
config/Phase_Porosity.config → examples/ConfigFiles/Phase_Porosity.config
View File

0
config/Phase_Thermal.config → examples/ConfigFiles/Phase_Thermal.config
View File

0
config/Phase_Vacancy.config → examples/ConfigFiles/Phase_Vacancy.config
View File

0
config/Source_Damage_IsoBrittle.config → examples/ConfigFiles/Source_Damage_IsoBrittle.config
View File

0
config/Source_Thermal_Dissipation.config → examples/ConfigFiles/Source_Thermal_Dissipation.config
View File

0
config/Source_Vacancy_Irradiation.config → examples/ConfigFiles/Source_Vacancy_Irradiation.config
View File

0
config/Source_Vacancy_PhenoPlasticity.config → examples/ConfigFiles/Source_Vacancy_PhenoPlasticity.config
View File

0
config/Texture_FiberExample.config → examples/ConfigFiles/Texture_FiberExample.config
View File

0
config/Texture_Gauss_001.config → examples/ConfigFiles/Texture_Gauss_001.config
View File

0
config/Texture_Gauss_101.config → examples/ConfigFiles/Texture_Gauss_101.config
View File

0
config/Texture_Gauss_111.config → examples/ConfigFiles/Texture_Gauss_111.config
View File

0
config/Texture_Gauss_123.config → examples/ConfigFiles/Texture_Gauss_123.config
View File

0
config/Texture_RandomSingleCrystals.config → examples/ConfigFiles/Texture_RandomSingleCrystals.config
View File

0
config/Texture_Rolling.config → examples/ConfigFiles/Texture_Rolling.config
View File

0
config/debug.config → examples/ConfigFiles/debug.config
View File

0
config/material.config → examples/ConfigFiles/material.config
View File

1
config/numerics.config → examples/ConfigFiles/numerics.config
View File

@ -8,7 +8,6 @@ pert_Fg 1.0e-7 # deformation gradient perturbation for g
pert_method 1 # perturbation method (1 = forward, 2 = backward or 3 = central) pert_method 1 # perturbation method (1 = forward, 2 = backward or 3 = central)
integrator 1 # integration method (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp) integrator 1 # integration method (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp)
integratorStiffness 1 # integration method used for stiffness (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp) integratorStiffness 1 # integration method used for stiffness (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp)
analyticJaco 1 # use analytic Jacobian or perturbation (0 = perturbations, 1 = analytic)
unitlength 1 # physical length of one computational length unit unitlength 1 # physical length of one computational length unit
usepingpong 1 # use the ping pong (collect <-> calc) scheme (always off for Abaqus exp, must be on for Spectral Solver) usepingpong 1 # use the ping pong (collect <-> calc) scheme (always off for Abaqus exp, must be on for Spectral Solver)

0
config/rollingTexture.linearODF → examples/ConfigFiles/rollingTexture.linearODF
View File

83
installation/mods_MarcMentat/apply_DAMASK_modifications.sh
View File

@ -1,21 +1,24 @@
#!/usr/bin/env bash #!/usr/bin/env bash
DEFAULT_VERSION='2015' SCRIPTLOCATION="$( cd "$( dirname "$0" )" && pwd )"
DAMASK_ROOT=$SCRIPTLOCATION/../../
WORKINGDIR="$( cd "$( dirname "$0" )" && pwd )" # defining set() allows to source the same file for tcsh and bash, with and without space around =
set() {
if [ -f $HOME/.damask/damask.conf ]; then export $1$2$3
source $HOME/.damask/damask.conf }
else source $DAMASK_ROOT/CONFIG
source /etc/damask.conf
fi
if [ "x$MSC_ROOT" != "x" ]; then if [ "x$MSC_ROOT" != "x" ]; then
DEFAULT_DIR=$MSC_ROOT DEFAULT_DIR=$MSC_ROOT
fi fi
if [ "x$MARC_VERSION" != "x" ]; then
DEFAULT_VERSION=$MARC_VERSION
fi
if [ "x$DAMASK_BIN" != "x" ]; then
BIN_DIR=$DAMASK_BIN
fi
while [ ! -d "$SCRIPTLOCATION/$VERSION" ] || [ -z "$VERSION" ]
while [ ! -d "$WORKINGDIR/$VERSION" ] || [ -z "$VERSION" ]
do do
echo "Input version of MARC/MENTAT installation: [${DEFAULT_VERSION}]" echo "Input version of MARC/MENTAT installation: [${DEFAULT_VERSION}]"
read VERSION read VERSION
@ -66,7 +69,7 @@ for filename in 'comp_damask' \
'run_damask_lmp' \ 'run_damask_lmp' \
'run_damask_hmp' \ 'run_damask_hmp' \
'include_linux64'; do 'include_linux64'; do
cp $WORKINGDIR/$VERSION/Marc_tools/$filename $theDIR cp $SCRIPTLOCATION/$VERSION/Marc_tools/$filename $theDIR
echo $theDIR/$filename | xargs perl -pi -e "s:%INSTALLDIR%:${INSTALLDIR}:g" echo $theDIR/$filename | xargs perl -pi -e "s:%INSTALLDIR%:${INSTALLDIR}:g"
echo $theDIR/$filename | xargs perl -pi -e "s:%VERSION%:${VERSION}:g" echo $theDIR/$filename | xargs perl -pi -e "s:%VERSION%:${VERSION}:g"
echo $filename echo $filename
@ -89,7 +92,7 @@ for filename in 'edit_window' \
'kill7' \ 'kill7' \
'kill8' \ 'kill8' \
'kill9'; do 'kill9'; do
cp $WORKINGDIR/$VERSION/Mentat_bin/$filename $theDIR cp $SCRIPTLOCATION/$VERSION/Mentat_bin/$filename $theDIR
echo $theDIR/$filename | xargs perl -pi -e "s:%INSTALLDIR%:${INSTALLDIR}:g" echo $theDIR/$filename | xargs perl -pi -e "s:%INSTALLDIR%:${INSTALLDIR}:g"
echo $theDIR/$filename | xargs perl -pi -e "s:%VERSION%:${VERSION}:g" echo $theDIR/$filename | xargs perl -pi -e "s:%VERSION%:${VERSION}:g"
echo $theDIR/$filename | xargs perl -pi -e "s:%EDITOR%:${EDITOR}:g" echo $theDIR/$filename | xargs perl -pi -e "s:%EDITOR%:${EDITOR}:g"
@ -101,7 +104,7 @@ echo ''
echo 'copying Mentat menus...' echo 'copying Mentat menus...'
theDIR=$INSTALLDIR/mentat$VERSION/menus theDIR=$INSTALLDIR/mentat$VERSION/menus
for filename in 'job_run.ms'; do for filename in 'job_run.ms'; do
cp $WORKINGDIR/$VERSION/Mentat_menus/$filename $theDIR cp $SCRIPTLOCATION/$VERSION/Mentat_menus/$filename $theDIR
echo $theDIR/$filename | xargs perl -pi -e "s:%INSTALLDIR%:${INSTALLDIR}:g" echo $theDIR/$filename | xargs perl -pi -e "s:%INSTALLDIR%:${INSTALLDIR}:g"
echo $theDIR/$filename | xargs perl -pi -e "s:%VERSION%:${VERSION}:g" echo $theDIR/$filename | xargs perl -pi -e "s:%VERSION%:${VERSION}:g"
echo $filename echo $filename
@ -121,31 +124,33 @@ chmod 755 $INSTALLDIR/mentat$VERSION/bin/submit{4..9}
chmod 755 $INSTALLDIR/mentat$VERSION/bin/kill{4..9} chmod 755 $INSTALLDIR/mentat$VERSION/bin/kill{4..9}
#creating symlinks for run_damask_scripts in /usr/local/bin #creating symlinks for run_damask_scripts in /usr/local/bin
BIN_DIR=/usr/local/bin
echo '' if [ -d "$BIN_DIR" ]; then
echo "Do you want to create symlinks for run_damask scripts in ${BIN_DIR} [YES/no] ?" echo ''
read YESNO echo "Do you want to create symlinks for run_damask scripts in ${BIN_DIR} [YES/no] ?"
if [ -z "$YESNO" ]; then read YESNO
YESNO=yes if [ -z "$YESNO" ]; then
fi YESNO=yes
case $YESNO in fi
y* | Y* ) case $YESNO in
echo'' y* | Y* )
echo 'creating symlinks ...' echo''
echo'' echo 'creating symlinks ...'
theDIR=$INSTALLDIR/marc$VERSION/tools echo''
for filename in 'run_damask' \ theDIR=$INSTALLDIR/marc$VERSION/tools
'run_damask_l' \ for filename in 'run_damask' \
'run_damask_h' \ 'run_damask_l' \
'run_damask_mp' \ 'run_damask_h' \
'run_damask_lmp' \ 'run_damask_mp' \
'run_damask_hmp'; do 'run_damask_lmp' \
echo ${filename:4}$VERSION 'run_damask_hmp'; do
[ -f $BIN_DIR/${filename:4}$VERSION ] && rm $BIN_DIR/${filename:4}$VERSION echo ${filename:4}$VERSION
ln -s $theDIR/$filename $BIN_DIR/${filename:4}$VERSION [ -f $BIN_DIR/${filename:4}$VERSION ] && rm $BIN_DIR/${filename:4}$VERSION
done ln -s $theDIR/$filename $BIN_DIR/${filename:4}$VERSION
;; done
esac ;;
esac
fi
echo '' echo ''
echo 'done.' echo 'done.'

14
installation/mods_MarcMentat/installation.txt
View File

@ -1,21 +1,13 @@
Install MPIE modifications to use DAMASK_marc2010/11/12/13/14 Install DAMASK modifications to use DAMASK_marc
This is for the Linux64 version of Marc/Mentat
Refer to http://damask.mpie.de for complete installation instructions. Refer to http://damask.mpie.de for complete installation instructions.
Usually you will need to be root for this to work! Usually you will need to be root for this to work!
This is for the Linux64 version of Marc/Mentat2011/12/13(.1)/14(.2).
See Marc and Mentat Release Guide for List of Build and Supported Platforms! See Marc and Mentat Release Guide for List of Build and Supported Platforms!
The Intel Fortran compiler needs to be installed. Marc/Mentat as well as DAMASK rely on Intel Fortran 12.0 or above! Make sure that ifort (the compiler executable) is in the PATH and that LD_LIBRARY_PATH is set correctly, refer to the Intel installation guide for instructions on how to do this. The Intel Fortran compiler needs to be installed.
The AMD Core Math Library or an other BLAS implementation (currently IMKL and LAPACK are supported) needs to be installed!
Add acml path to LD_LIBRARY_PATH, to do so either use the script setup_shellrc.py in the installation folder or for a system wide setup edit /etc/csh.cshrc.local and/or /etc/bash.bashrc.local.
Assuming ACML is installed in path ACMLDIR the path should read: /ACMLDIR/ifort64_mp/lib (this is the version using paralllization, ie openmp):
for bash: LD_LIBRARY_PATH="/ACMLDIR/ifort64_mp/lib:${LD_LIBRARY_PATH}"; export LD_LIBRARY_PATH
for csh: setenv LD_LIBRARY_PATH /ACMLDIR/ifort64_mp/lib:${LD_LIBRARY_PATH}
As the apply_DAMASK_modifications script has to fix the path for the BLAS it needs to be installed in a place that can be accessed by all users of the system. In addition you have to specify the respective locations in DAMASK_ROOT/lib/pathInfo.
1) Install Marc, Mentat and Documentation as usual 1) Install Marc, Mentat and Documentation as usual
Run the test example including subroutines to confirm that the installation of both Marc/Mentat and the Intel Fortran Compiler is ok! Run the test example including subroutines to confirm that the installation of both Marc/Mentat and the Intel Fortran Compiler is ok!

73
installation/symlink_Code.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2 #!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys import os,sys
@ -11,39 +11,58 @@ bin_link = { \
], ],
} }
MarcReleases =[2011,2012,2013,2013.1,2014,2014.2,2015] MarcReleases =[ \
'2011',
'2012',
'2013',
'2013.1',
'2014',
'2014.2',
'2015',
]
baseDir = damask.Environment('../../').relPath('code/') damaskEnv = damask.Environment()
baseDir = damaskEnv.relPath('code/')
try: binDir = damaskEnv.options['DAMASK_BIN']
binDir = damask.Environment().options['DAMASK_BIN']
except:
root=os.access('/usr/local/bin', os.W_OK)
if root:
binDir = '/usr/local/bin'
else:
binDir = os.path.join(os.getenv('HOME'),'bin')
if not os.path.isdir(binDir): if not os.path.isdir(binDir):
os.mkdir(binDir) os.mkdir(binDir)
for dir in bin_link: sys.stdout.write('\nsymbolic linking...\n')
for file in bin_link[dir]:
src = os.path.abspath(os.path.join(baseDir,dir,file))
if os.path.exists(src):
sym_link = os.path.abspath(os.path.join(binDir,\
{True: dir,
False:os.path.splitext(file)[0]}[file == '']))
if os.path.lexists(sym_link): os.remove(sym_link)
os.symlink(src,sym_link)
sys.stdout.write(sym_link+' -> '+src+'\n')
for subDir in bin_link:
theDir = os.path.abspath(os.path.join(baseDir,subDir))
sys.stdout.write('\n'+binDir+' ->\n'+theDir+damask.util.deemph(' ...')+'\n')
for theFile in bin_link[subDir]:
theName,theExt = os.path.splitext(theFile)
src = os.path.abspath(os.path.join(theDir,theFile))
if os.path.exists(src):
sym_link = os.path.abspath(os.path.join(binDir,subDir if theFile == '' else theName))
if os.path.lexists(sym_link):
os.remove(sym_link)
output = theName+damask.util.deemph(theExt)
else:
output = damask.util.emph(theName)+damask.util.deemph(theExt)
sys.stdout.write(damask.util.deemph('... ')+output+'\n')
os.symlink(src,sym_link)
sys.stdout.write('\nMSC.Marc versioning...\n\n')
theMaster = 'DAMASK_marc.f90'
for version in MarcReleases: for version in MarcReleases:
src = os.path.abspath(os.path.join(baseDir,'DAMASK_marc.f90')) src = os.path.abspath(os.path.join(baseDir,theMaster))
if os.path.exists(src): if os.path.exists(src):
sym_link = os.path.abspath(os.path.join(baseDir,'DAMASK_marc'+str(version)+'.f90')) sym_link = os.path.abspath(os.path.join(baseDir,'DAMASK_marc{}.f90'.format(version)))
if os.path.lexists(sym_link): os.remove(sym_link) if os.path.lexists(sym_link):
os.symlink(os.path.relpath(src,baseDir),sym_link) os.remove(sym_link)
sys.stdout.write(sym_link+' -> '+src+'\n') output = version
else:
output = damask.util.emph(version)
sys.stdout.write(' '+output+'\n')
os.symlink(theMaster,sym_link)

Some files were not shown because too many files have changed in this diff Show More