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:
Philip Eisenlohr 2018-08-31 11:23:25 -04:00
parent 3e81110830 b24ebb8a5b
commit 6a64637ac3
55 changed files with 3369 additions and 26907 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
.gitlab-ci.yml
View File

@ -3,8 +3,8 @@ stages:
- prepareAll - prepareAll
- preprocessing - preprocessing
- postprocessing - postprocessing
- compileSpectralIntel - compilePETScIntel
- compileSpectralGNU - compilePETScGNU
- prepareSpectral - prepareSpectral
- spectral - spectral
- compileMarc2017 - compileMarc2017
@ -186,8 +186,8 @@ Post_ParaviewRelated:
- release - release
################################################################################################### ###################################################################################################
Compile_Intel: Compile_Spectral_Intel:
stage: compileSpectralIntel stage: compilePETScIntel
script: script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- SpectralAll_compile/test.py - SpectralAll_compile/test.py
@ -195,9 +195,18 @@ Compile_Intel:
- master - master
- release - release
Compile_FEM_Intel:
stage: compilePETScIntel
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- FEM_compile/test.py
except:
- master
- release
################################################################################################### ###################################################################################################
Compile_GNU: Compile_Spectral_GNU:
stage: compileSpectralGNU stage: compilePETScGNU
script: script:
- module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU - module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU
- SpectralAll_compile/test.py - SpectralAll_compile/test.py
@ -205,6 +214,15 @@ Compile_GNU:
- master - master
- release - release
Compile_FEM_GNU:
stage: compilePETScGNU
script:
- module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU
- FEM_compile/test.py
except:
- master
- release
################################################################################################### ###################################################################################################
Compile_Intel_Prepare: Compile_Intel_Prepare:
stage: prepareSpectral stage: prepareSpectral

2
PRIVATE

@ -1 +1 @@
Subproject commit c4471725893e301044924eb0990e2ad619aa0a46 Subproject commit fa02113fa7a0af3376648e4320318ec337fe79aa

2
VERSION
View File

@ -1 +1 @@
v2.0.2-381-gc03ea8f5 v2.0.2-490-g29e55d20

2
examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Martensite.config
View File

@ -18,5 +18,5 @@ tau0_slip 405.8e6 456.7e6 # per family
tausat_slip 872.9e6 971.2e6 # per family tausat_slip 872.9e6 971.2e6 # per family
h0_slipslip 563.0e9 h0_slipslip 563.0e9
interaction_slipslip 1 1 1.4 1.4 1.4 1.4 interaction_slipslip 1 1 1.4 1.4 1.4 1.4
w0_slip 2.0 a_slip 2.0
(output) totalshear (output) totalshear

52
installation/mods_MarcMentat/2016/Marc_tools/comp_damask
View File

@ -1,52 +0,0 @@
#!/bin/ksh
# 1st arg: $DIR
# 2nd arg: $DIRJOB
# 3rd arg: $user
# 4th arg: $program
DIR=$1
user=$3
program=$4
usernoext=$user
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
# add BLAS options for linking
BLAS="%BLAS%"
. $DIR/tools/include
DIRJOB=$2
cd $DIRJOB
echo "Compiling and linking user subroutine $user on host `hostname`"
echo "program: $program"
$DFORTRAN $user || \
{
echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
userobj=$usernoext.o
$LOAD ${program} $DIR/lib/main.o\
$DIR/lib/blkdta.o $DIR/lib/comm?.o \
${userobj-} \
$DIR/lib/srclib.a \
$MNFLIBS \
$MDUSER \
../lib/mdsrc.a \
../lib/mcvfit.a \
$STUBS \
${SOLVERLIBS} \
$TKLIBS \
$MRCLIBS \
$METISLIBS \
$BLAS \
$SYSLIBS || \
{
echo "$0: link failed for $usernoext.o on host `hostname`"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod

52
installation/mods_MarcMentat/2016/Marc_tools/comp_damask_h
View File

@ -1,52 +0,0 @@
#!/bin/ksh
# 1st arg: $DIR
# 2nd arg: $DIRJOB
# 3rd arg: $user
# 4th arg: $program
DIR=$1
user=$3
program=$4
usernoext=$user
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
# add BLAS options for linking
BLAS="%BLAS%"
. $DIR/tools/include
DIRJOB=$2
cd $DIRJOB
echo "Compiling and linking user subroutine $user on host `hostname`"
echo "program: $program"
$DFORTHIGH $user || \
{
echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
userobj=$usernoext.o
$LOAD ${program} $DIR/lib/main.o\
$DIR/lib/blkdta.o $DIR/lib/comm?.o \
${userobj-} \
$DIR/lib/srclib.a \
$MNFLIBS \
$MDUSER \
../lib/mdsrc.a \
../lib/mcvfit.a \
$STUBS \
${SOLVERLIBS} \
$TKLIBS \
$MRCLIBS \
$METISLIBS \
$BLAS \
$SYSLIBS || \
{
echo "$0: link failed for $usernoext.o on host `hostname`"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod

52
installation/mods_MarcMentat/2016/Marc_tools/comp_damask_l
View File

@ -1,52 +0,0 @@
#!/bin/ksh
# 1st arg: $DIR
# 2nd arg: $DIRJOB
# 3rd arg: $user
# 4th arg: $program
DIR=$1
user=$3
program=$4
usernoext=$user
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
# add BLAS options for linking
BLAS="%BLAS%"
. $DIR/tools/include
DIRJOB=$2
cd $DIRJOB
echo "Compiling and linking user subroutine $user on host `hostname`"
echo "program: $program"
$DFORTLOW $user || \
{
echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
userobj=$usernoext.o
$LOAD ${program} $DIR/lib/main.o\
$DIR/lib/blkdta.o $DIR/lib/comm?.o \
${userobj-} \
$DIR/lib/srclib.a \
$MNFLIBS \
$MDUSER \
../lib/mdsrc.a \
../lib/mcvfit.a \
$STUBS \
${SOLVERLIBS} \
$TKLIBS \
$MRCLIBS \
$METISLIBS \
$BLAS \
$SYSLIBS || \
{
echo "$0: link failed for $usernoext.o on host `hostname`"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod

4112
installation/mods_MarcMentat/2016/Marc_tools/run_damask
View File

File diff suppressed because it is too large Load Diff

4112
installation/mods_MarcMentat/2016/Marc_tools/run_damask_h
View File

File diff suppressed because it is too large Load Diff

4112
installation/mods_MarcMentat/2016/Marc_tools/run_damask_l
View File

File diff suppressed because it is too large Load Diff

8
installation/mods_MarcMentat/2016/Mentat_bin/kill7
View File

@ -1,8 +0,0 @@
#!/bin/sh
if [ "$1" = "" ]; then
echo "usage: $0 job_name"
exit 1
fi
echo STOP > $1.cnt

8
installation/mods_MarcMentat/2016/Mentat_bin/kill8
View File

@ -1,8 +0,0 @@
#!/bin/sh
if [ "$1" = "" ]; then
echo "usage: $0 job_name"
exit 1
fi
echo STOP > $1.cnt

8
installation/mods_MarcMentat/2016/Mentat_bin/kill9
View File

@ -1,8 +0,0 @@
#!/bin/sh
if [ "$1" = "" ]; then
echo "usage: $0 job_name"
exit 1
fi
echo STOP > $1.cnt

187
installation/mods_MarcMentat/2016/Mentat_bin/submit7
View File

@ -1,187 +0,0 @@
#!/bin/sh
#
# The exit status of this script is read by Mentat.
# Normal exit status is 0.
#
DIR=%INSTALLDIR%/marc%VERSION%
if test $MARCDIR1
then
DIR=$MARCDIR1
fi
if test -z "$DIR"; then
REALCOM="`ls -l $0 |awk '{ print $NF; }'`"
DIRSCRIPT=`dirname $REALCOM`
case $DIRSCRIPT in
\/*)
;;
*)
DIRSCRIPT=`pwd`/$DIRSCRIPT
;;
esac
. $DIRSCRIPT/getarch
DIR="$MENTAT_MARCDIR"
fi
SRCEXT=.f
SRCEXTC=.F
RSTEXT=.t08
PSTEXT=.t19
PSTEXTB=.t16
VWFCEXT=.vfs
slv=$1
version=$2
ndom_fea_solver=$3
ndom_preprocessor=$4
hostfile=$5
compat=$6
job=$7
srcfile=$8
srcmeth=$9
shift 9 # cannot use $10, $11, ...
restart=$1
postfile=$2
viewfactorsfile=$3
autorst=$4
copy_datfile="-ci $5"
copy_postfile="-cr $6"
scr_dir=$7
dcoup=$8
assem_recov_nthread=$9
shift 9 # cannot use $10, $11, ...
nthread=$1
nsolver=$2
mode=$3
gpu=$4
if [ "$slv" != "" -a "$slv" != "marc" ]; then
slv="-iam sfm"
else
slv=""
fi
if [ "$ndom_fea_solver" != "" -a "$ndom_fea_solver" != "1" ]; then
nprocds="-nprocds $ndom_fea_solver"
else
nprocd=""
if [ "$ndom_preprocessor" != "" -a "$ndom_preprocessor" != "1" ]; then
nprocd="-nprocd $ndom_preprocessor"
else
nprocd=""
fi
fi
if [ "$srcfile" != "" -a "$srcfile" != "-" ]; then
srcfile=`echo $srcfile | sed "s/$SRCEXT$//" | sed "s/$SRCEXTC$//"`
case "$srcmeth" in
-)
srcfile="-u $srcfile"
;;
compsave)
srcfile="-u $srcfile -save y"
;;
runsaved)
srcfile=${srcfile%.*}".marc"
srcfile="-prog $srcfile"
;;
esac
else
srcfile=""
fi
if [ "$restart" != "" -a "$restart" != "-" ]; then
restart=`echo $restart | sed "s/$RSTEXT$//"`
restart="-r $restart"
else
restart=""
fi
if [ "$postfile" != "" -a "$postfile" != "-" ]; then
postfile=`echo $postfile | sed "s/$PSTEXT$//"`
postfile=`echo $postfile | sed "s/$PSTEXTB$//"`
postfile="-pid $postfile"
else
postfile=""
fi
if [ "$viewfactorsfile" != "" -a "$viewfactorsfile" != "-" ]; then
viewfactorsfile=`echo $viewfactorsfile | sed "s/$VWFCEXT$//"`
viewfactorsfile="-vf $viewfactorsfile"
else
viewfactorsfile=""
fi
if [ "$hostfile" != "" -a "$hostfile" != "-" ]; then
hostfile="-ho $hostfile"
else
hostfile=""
fi
if [ "$compat" != "" -a "$compat" != "-" ]; then
compat="-co $compat"
else
compat=""
fi
if [ "$scr_dir" != "" -a "$scr_dir" != "-" ]; then
scr_dir="-sd $scr_dir"
else
scr_dir=""
fi
if [ "$dcoup" != "" -a "$dcoup" != "0" ]; then
dcoup="-dcoup $dcoup"
else
dcoup=""
fi
if [ "$assem_recov_nthread" != "" -a "$assem_recov_nthread" != "1" ]; then
assem_recov_nthread="-nthread_elem $assem_recov_nthread"
else
assem_recov_nthread=""
fi
if [ "$nthread" != "" -a "$nthread" != "0" -a "$nthread" != "1" ]; then
nthread="-nthread $nthread"
else
nthread=""
fi
if [ "$nsolver" != "" -a "$nsolver" != "0" ]; then
nsolver="-nsolver $nsolver"
else
nsolver=""
fi
case "$mode" in
4) mode="-mo i4" ;;
8) mode="-mo i8" ;;
*) mode= ;;
esac
if [ "$gpu" != "" -a "$gpu" != "-" ]; then
gpu="-gpu $gpu"
else
gpu=""
fi
rm -f $job.cnt
rm -f $job.sts
rm -f $job.out
rm -f $job.log
# To prevent a mismatch with the python version used by the solver
# do *not* prepend $MENTAT_INSTALL_DIR/python/bin to environment variable PATH
# unset environment variables PYTHONHOME and PYTHONPATH
unset PYTHONHOME
unset PYTHONPATH
"${DIR}/tools/run_damask_h" $slv -j $job -v n -b y $nprocds $nprocd -autorst $autorst \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu > /dev/null 2>&1
sleep 1
exit 0

187
installation/mods_MarcMentat/2016/Mentat_bin/submit8
View File

@ -1,187 +0,0 @@
#!/bin/sh
#
# The exit status of this script is read by Mentat.
# Normal exit status is 0.
#
DIR=%INSTALLDIR%/marc%VERSION%
if test $MARCDIR1
then
DIR=$MARCDIR1
fi
if test -z "$DIR"; then
REALCOM="`ls -l $0 |awk '{ print $NF; }'`"
DIRSCRIPT=`dirname $REALCOM`
case $DIRSCRIPT in
\/*)
;;
*)
DIRSCRIPT=`pwd`/$DIRSCRIPT
;;
esac
. $DIRSCRIPT/getarch
DIR="$MENTAT_MARCDIR"
fi
SRCEXT=.f
SRCEXTC=.F
RSTEXT=.t08
PSTEXT=.t19
PSTEXTB=.t16
VWFCEXT=.vfs
slv=$1
version=$2
ndom_fea_solver=$3
ndom_preprocessor=$4
hostfile=$5
compat=$6
job=$7
srcfile=$8
srcmeth=$9
shift 9 # cannot use $10, $11, ...
restart=$1
postfile=$2
viewfactorsfile=$3
autorst=$4
copy_datfile="-ci $5"
copy_postfile="-cr $6"
scr_dir=$7
dcoup=$8
assem_recov_nthread=$9
shift 9 # cannot use $10, $11, ...
nthread=$1
nsolver=$2
mode=$3
gpu=$4
if [ "$slv" != "" -a "$slv" != "marc" ]; then
slv="-iam sfm"
else
slv=""
fi
if [ "$ndom_fea_solver" != "" -a "$ndom_fea_solver" != "1" ]; then
nprocds="-nprocds $ndom_fea_solver"
else
nprocd=""
if [ "$ndom_preprocessor" != "" -a "$ndom_preprocessor" != "1" ]; then
nprocd="-nprocd $ndom_preprocessor"
else
nprocd=""
fi
fi
if [ "$srcfile" != "" -a "$srcfile" != "-" ]; then
srcfile=`echo $srcfile | sed "s/$SRCEXT$//" | sed "s/$SRCEXTC$//"`
case "$srcmeth" in
-)
srcfile="-u $srcfile"
;;
compsave)
srcfile="-u $srcfile -save y"
;;
runsaved)
srcfile=${srcfile%.*}".marc"
srcfile="-prog $srcfile"
;;
esac
else
srcfile=""
fi
if [ "$restart" != "" -a "$restart" != "-" ]; then
restart=`echo $restart | sed "s/$RSTEXT$//"`
restart="-r $restart"
else
restart=""
fi
if [ "$postfile" != "" -a "$postfile" != "-" ]; then
postfile=`echo $postfile | sed "s/$PSTEXT$//"`
postfile=`echo $postfile | sed "s/$PSTEXTB$//"`
postfile="-pid $postfile"
else
postfile=""
fi
if [ "$viewfactorsfile" != "" -a "$viewfactorsfile" != "-" ]; then
viewfactorsfile=`echo $viewfactorsfile | sed "s/$VWFCEXT$//"`
viewfactorsfile="-vf $viewfactorsfile"
else
viewfactorsfile=""
fi
if [ "$hostfile" != "" -a "$hostfile" != "-" ]; then
hostfile="-ho $hostfile"
else
hostfile=""
fi
if [ "$compat" != "" -a "$compat" != "-" ]; then
compat="-co $compat"
else
compat=""
fi
if [ "$scr_dir" != "" -a "$scr_dir" != "-" ]; then
scr_dir="-sd $scr_dir"
else
scr_dir=""
fi
if [ "$dcoup" != "" -a "$dcoup" != "0" ]; then
dcoup="-dcoup $dcoup"
else
dcoup=""
fi
if [ "$assem_recov_nthread" != "" -a "$assem_recov_nthread" != "1" ]; then
assem_recov_nthread="-nthread_elem $assem_recov_nthread"
else
assem_recov_nthread=""
fi
if [ "$nthread" != "" -a "$nthread" != "0" -a "$nthread" != "1" ]; then
nthread="-nthread $nthread"
else
nthread=""
fi
if [ "$nsolver" != "" -a "$nsolver" != "0" ]; then
nsolver="-nsolver $nsolver"
else
nsolver=""
fi
case "$mode" in
4) mode="-mo i4" ;;
8) mode="-mo i8" ;;
*) mode= ;;
esac
if [ "$gpu" != "" -a "$gpu" != "-" ]; then
gpu="-gpu $gpu"
else
gpu=""
fi
rm -f $job.cnt
rm -f $job.sts
rm -f $job.out
rm -f $job.log
# To prevent a mismatch with the python version used by the solver
# do *not* prepend $MENTAT_INSTALL_DIR/python/bin to environment variable PATH
# unset environment variables PYTHONHOME and PYTHONPATH
unset PYTHONHOME
unset PYTHONPATH
"${DIR}/tools/run_damask" $slv -j $job -v n -b y $nprocds $nprocd -autorst $autorst \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu > /dev/null 2>&1
sleep 1
exit 0

187
installation/mods_MarcMentat/2016/Mentat_bin/submit9
View File

@ -1,187 +0,0 @@
#!/bin/sh
#
# The exit status of this script is read by Mentat.
# Normal exit status is 0.
#
DIR=%INSTALLDIR%/marc%VERSION%
if test $MARCDIR1
then
DIR=$MARCDIR1
fi
if test -z "$DIR"; then
REALCOM="`ls -l $0 |awk '{ print $NF; }'`"
DIRSCRIPT=`dirname $REALCOM`
case $DIRSCRIPT in
\/*)
;;
*)
DIRSCRIPT=`pwd`/$DIRSCRIPT
;;
esac
. $DIRSCRIPT/getarch
DIR="$MENTAT_MARCDIR"
fi
SRCEXT=.f
SRCEXTC=.F
RSTEXT=.t08
PSTEXT=.t19
PSTEXTB=.t16
VWFCEXT=.vfs
slv=$1
version=$2
ndom_fea_solver=$3
ndom_preprocessor=$4
hostfile=$5
compat=$6
job=$7
srcfile=$8
srcmeth=$9
shift 9 # cannot use $10, $11, ...
restart=$1
postfile=$2
viewfactorsfile=$3
autorst=$4
copy_datfile="-ci $5"
copy_postfile="-cr $6"
scr_dir=$7
dcoup=$8
assem_recov_nthread=$9
shift 9 # cannot use $10, $11, ...
nthread=$1
nsolver=$2
mode=$3
gpu=$4
if [ "$slv" != "" -a "$slv" != "marc" ]; then
slv="-iam sfm"
else
slv=""
fi
if [ "$ndom_fea_solver" != "" -a "$ndom_fea_solver" != "1" ]; then
nprocds="-nprocds $ndom_fea_solver"
else
nprocd=""
if [ "$ndom_preprocessor" != "" -a "$ndom_preprocessor" != "1" ]; then
nprocd="-nprocd $ndom_preprocessor"
else
nprocd=""
fi
fi
if [ "$srcfile" != "" -a "$srcfile" != "-" ]; then
srcfile=`echo $srcfile | sed "s/$SRCEXT$//" | sed "s/$SRCEXTC$//"`
case "$srcmeth" in
-)
srcfile="-u $srcfile"
;;
compsave)
srcfile="-u $srcfile -save y"
;;
runsaved)
srcfile=${srcfile%.*}".marc"
srcfile="-prog $srcfile"
;;
esac
else
srcfile=""
fi
if [ "$restart" != "" -a "$restart" != "-" ]; then
restart=`echo $restart | sed "s/$RSTEXT$//"`
restart="-r $restart"
else
restart=""
fi
if [ "$postfile" != "" -a "$postfile" != "-" ]; then
postfile=`echo $postfile | sed "s/$PSTEXT$//"`
postfile=`echo $postfile | sed "s/$PSTEXTB$//"`
postfile="-pid $postfile"
else
postfile=""
fi
if [ "$viewfactorsfile" != "" -a "$viewfactorsfile" != "-" ]; then
viewfactorsfile=`echo $viewfactorsfile | sed "s/$VWFCEXT$//"`
viewfactorsfile="-vf $viewfactorsfile"
else
viewfactorsfile=""
fi
if [ "$hostfile" != "" -a "$hostfile" != "-" ]; then
hostfile="-ho $hostfile"
else
hostfile=""
fi
if [ "$compat" != "" -a "$compat" != "-" ]; then
compat="-co $compat"
else
compat=""
fi
if [ "$scr_dir" != "" -a "$scr_dir" != "-" ]; then
scr_dir="-sd $scr_dir"
else
scr_dir=""
fi
if [ "$dcoup" != "" -a "$dcoup" != "0" ]; then
dcoup="-dcoup $dcoup"
else
dcoup=""
fi
if [ "$assem_recov_nthread" != "" -a "$assem_recov_nthread" != "1" ]; then
assem_recov_nthread="-nthread_elem $assem_recov_nthread"
else
assem_recov_nthread=""
fi
if [ "$nthread" != "" -a "$nthread" != "0" -a "$nthread" != "1" ]; then
nthread="-nthread $nthread"
else
nthread=""
fi
if [ "$nsolver" != "" -a "$nsolver" != "0" ]; then
nsolver="-nsolver $nsolver"
else
nsolver=""
fi
case "$mode" in
4) mode="-mo i4" ;;
8) mode="-mo i8" ;;
*) mode= ;;
esac
if [ "$gpu" != "" -a "$gpu" != "-" ]; then
gpu="-gpu $gpu"
else
gpu=""
fi
rm -f $job.cnt
rm -f $job.sts
rm -f $job.out
rm -f $job.log
# To prevent a mismatch with the python version used by the solver
# do *not* prepend $MENTAT_INSTALL_DIR/python/bin to environment variable PATH
# unset environment variables PYTHONHOME and PYTHONPATH
unset PYTHONHOME
unset PYTHONPATH
"${DIR}/tools/run_damask_l" $slv -j $job -v n -b y $nprocds $nprocd -autorst $autorst \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu > /dev/null 2>&1
sleep 1
exit 0

52
installation/mods_MarcMentat/2017/Marc_tools/comp_damask
View File

@ -1,52 +0,0 @@
#!/bin/ksh
# 1st arg: $DIR
# 2nd arg: $DIRJOB
# 3rd arg: $user
# 4th arg: $program
DIR=$1
user=$3
program=$4
usernoext=$user
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
# add BLAS options for linking
BLAS="%BLAS%"
. $DIR/tools/include
DIRJOB=$2
cd $DIRJOB
echo "Compiling and linking user subroutine $user on host `hostname`"
echo "program: $program"
$DFORTRAN $user || \
{
echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
userobj=$usernoext.o
$LOAD ${program} $DIR/lib/main.o\
$DIR/lib/blkdta.o $DIR/lib/comm?.o \
${userobj-} \
$DIR/lib/srclib.a \
$MNFLIBS \
$MDUSER \
../lib/mdsrc.a \
../lib/mcvfit.a \
$STUBS \
${SOLVERLIBS} \
$TKLIBS \
$MRCLIBS \
$METISLIBS \
$BLAS \
$SYSLIBS || \
{
echo "$0: link failed for $usernoext.o on host `hostname`"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod

52
installation/mods_MarcMentat/2017/Marc_tools/comp_damask_h
View File

@ -1,52 +0,0 @@
#!/bin/ksh
# 1st arg: $DIR
# 2nd arg: $DIRJOB
# 3rd arg: $user
# 4th arg: $program
DIR=$1
user=$3
program=$4
usernoext=$user
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
# add BLAS options for linking
BLAS="%BLAS%"
. $DIR/tools/include
DIRJOB=$2
cd $DIRJOB
echo "Compiling and linking user subroutine $user on host `hostname`"
echo "program: $program"
$DFORTHIGH $user || \
{
echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
userobj=$usernoext.o
$LOAD ${program} $DIR/lib/main.o\
$DIR/lib/blkdta.o $DIR/lib/comm?.o \
${userobj-} \
$DIR/lib/srclib.a \
$MNFLIBS \
$MDUSER \
../lib/mdsrc.a \
../lib/mcvfit.a \
$STUBS \
${SOLVERLIBS} \
$TKLIBS \
$MRCLIBS \
$METISLIBS \
$BLAS \
$SYSLIBS || \
{
echo "$0: link failed for $usernoext.o on host `hostname`"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod

52
installation/mods_MarcMentat/2017/Marc_tools/comp_damask_l
View File

@ -1,52 +0,0 @@
#!/bin/ksh
# 1st arg: $DIR
# 2nd arg: $DIRJOB
# 3rd arg: $user
# 4th arg: $program
DIR=$1
user=$3
program=$4
usernoext=$user
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .for`
usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f90`
# add BLAS options for linking
BLAS="%BLAS%"
. $DIR/tools/include
DIRJOB=$2
cd $DIRJOB
echo "Compiling and linking user subroutine $user on host `hostname`"
echo "program: $program"
$DFORTLOW $user || \
{
echo "$0: compile failed for $user"
exit 1
}
/bin/rm $program 2>/dev/null
userobj=$usernoext.o
$LOAD ${program} $DIR/lib/main.o\
$DIR/lib/blkdta.o $DIR/lib/comm?.o \
${userobj-} \
$DIR/lib/srclib.a \
$MNFLIBS \
$MDUSER \
../lib/mdsrc.a \
../lib/mcvfit.a \
$STUBS \
${SOLVERLIBS} \
$TKLIBS \
$MRCLIBS \
$METISLIBS \
$BLAS \
$SYSLIBS || \
{
echo "$0: link failed for $usernoext.o on host `hostname`"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod

4122
installation/mods_MarcMentat/2017/Marc_tools/run_damask
View File

File diff suppressed because it is too large Load Diff

4122
installation/mods_MarcMentat/2017/Marc_tools/run_damask_h
View File

File diff suppressed because it is too large Load Diff

4122
installation/mods_MarcMentat/2017/Marc_tools/run_damask_l
View File

File diff suppressed because it is too large Load Diff

8
installation/mods_MarcMentat/2017/Mentat_bin/kill7
View File

@ -1,8 +0,0 @@
#!/bin/sh
if [ "$1" = "" ]; then
echo "usage: $0 job_name"
exit 1
fi
echo STOP > $1.cnt

8
installation/mods_MarcMentat/2017/Mentat_bin/kill8
View File

@ -1,8 +0,0 @@
#!/bin/sh
if [ "$1" = "" ]; then
echo "usage: $0 job_name"
exit 1
fi
echo STOP > $1.cnt

8
installation/mods_MarcMentat/2017/Mentat_bin/kill9
View File

@ -1,8 +0,0 @@
#!/bin/sh
if [ "$1" = "" ]; then
echo "usage: $0 job_name"
exit 1
fi
echo STOP > $1.cnt

187
installation/mods_MarcMentat/2017/Mentat_bin/submit7
View File

@ -1,187 +0,0 @@
#!/bin/sh
#
# The exit status of this script is read by Mentat.
# Normal exit status is 0.
#
DIR=%INSTALLDIR%/marc%VERSION%
if test $MARCDIR1
then
DIR=$MARCDIR1
fi
if test -z "$DIR"; then
REALCOM="`ls -l $0 |awk '{ print $NF; }'`"
DIRSCRIPT=`dirname $REALCOM`
case $DIRSCRIPT in
\/*)
;;
*)
DIRSCRIPT=`pwd`/$DIRSCRIPT
;;
esac
. $DIRSCRIPT/getarch
DIR="$MENTAT_MARCDIR"
fi
SRCEXT=.f
SRCEXTC=.F
RSTEXT=.t08
PSTEXT=.t19
PSTEXTB=.t16
VWFCEXT=.vfs
slv=$1
version=$2
ndom_fea_solver=$3
ndom_preprocessor=$4
hostfile=$5
compat=$6
job=$7
srcfile=$8
srcmeth=$9
shift 9 # cannot use $10, $11, ...
restart=$1
postfile=$2
viewfactorsfile=$3
autorst=$4
copy_datfile="-ci $5"
copy_postfile="-cr $6"
scr_dir=$7
dcoup=$8
assem_recov_nthread=$9
shift 9 # cannot use $10, $11, ...
nthread=$1
nsolver=$2
mode=$3
gpu=$4
if [ "$slv" != "" -a "$slv" != "marc" ]; then
slv="-iam sfm"
else
slv=""
fi
if [ "$ndom_fea_solver" != "" -a "$ndom_fea_solver" != "1" ]; then
nprocds="-nprocds $ndom_fea_solver"
else
nprocd=""
if [ "$ndom_preprocessor" != "" -a "$ndom_preprocessor" != "1" ]; then
nprocd="-nprocd $ndom_preprocessor"
else
nprocd=""
fi
fi
if [ "$srcfile" != "" -a "$srcfile" != "-" ]; then
srcfile=`echo $srcfile | sed "s/$SRCEXT$//" | sed "s/$SRCEXTC$//"`
case "$srcmeth" in
-)
srcfile="-u $srcfile"
;;
compsave)
srcfile="-u $srcfile -save y"
;;
runsaved)
srcfile=${srcfile%.*}".marc"
srcfile="-prog $srcfile"
;;
esac
else
srcfile=""
fi
if [ "$restart" != "" -a "$restart" != "-" ]; then
restart=`echo $restart | sed "s/$RSTEXT$//"`
restart="-r $restart"
else
restart=""
fi
if [ "$postfile" != "" -a "$postfile" != "-" ]; then
postfile=`echo $postfile | sed "s/$PSTEXT$//"`
postfile=`echo $postfile | sed "s/$PSTEXTB$//"`
postfile="-pid $postfile"
else
postfile=""
fi
if [ "$viewfactorsfile" != "" -a "$viewfactorsfile" != "-" ]; then
viewfactorsfile=`echo $viewfactorsfile | sed "s/$VWFCEXT$//"`
viewfactorsfile="-vf $viewfactorsfile"
else
viewfactorsfile=""
fi
if [ "$hostfile" != "" -a "$hostfile" != "-" ]; then
hostfile="-ho $hostfile"
else
hostfile=""
fi
if [ "$compat" != "" -a "$compat" != "-" ]; then
compat="-co $compat"
else
compat=""
fi
if [ "$scr_dir" != "" -a "$scr_dir" != "-" ]; then
scr_dir="-sd $scr_dir"
else
scr_dir=""
fi
if [ "$dcoup" != "" -a "$dcoup" != "0" ]; then
dcoup="-dcoup $dcoup"
else
dcoup=""
fi
if [ "$assem_recov_nthread" != "" -a "$assem_recov_nthread" != "1" ]; then
assem_recov_nthread="-nthread_elem $assem_recov_nthread"
else
assem_recov_nthread=""
fi
if [ "$nthread" != "" -a "$nthread" != "0" -a "$nthread" != "1" ]; then
nthread="-nthread $nthread"
else
nthread=""
fi
if [ "$nsolver" != "" -a "$nsolver" != "0" ]; then
nsolver="-nsolver $nsolver"
else
nsolver=""
fi
case "$mode" in
4) mode="-mo i4" ;;
8) mode="-mo i8" ;;
*) mode= ;;
esac
if [ "$gpu" != "" -a "$gpu" != "-" ]; then
gpu="-gpu $gpu"
else
gpu=""
fi
rm -f $job.cnt
rm -f $job.sts
rm -f $job.out
rm -f $job.log
# To prevent a mismatch with the python version used by the solver
# do *not* prepend $MENTAT_INSTALL_DIR/python/bin to environment variable PATH
# unset environment variables PYTHONHOME and PYTHONPATH
unset PYTHONHOME
unset PYTHONPATH
"${DIR}/tools/run_damask_h" $slv -j $job -v n -b y $nprocds $nprocd -autorst $autorst \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu > /dev/null 2>&1
sleep 1
exit 0

187
installation/mods_MarcMentat/2017/Mentat_bin/submit8
View File

@ -1,187 +0,0 @@
#!/bin/sh
#
# The exit status of this script is read by Mentat.
# Normal exit status is 0.
#
DIR=%INSTALLDIR%/marc%VERSION%
if test $MARCDIR1
then
DIR=$MARCDIR1
fi
if test -z "$DIR"; then
REALCOM="`ls -l $0 |awk '{ print $NF; }'`"
DIRSCRIPT=`dirname $REALCOM`
case $DIRSCRIPT in
\/*)
;;
*)
DIRSCRIPT=`pwd`/$DIRSCRIPT
;;
esac
. $DIRSCRIPT/getarch
DIR="$MENTAT_MARCDIR"
fi
SRCEXT=.f
SRCEXTC=.F
RSTEXT=.t08
PSTEXT=.t19
PSTEXTB=.t16
VWFCEXT=.vfs
slv=$1
version=$2
ndom_fea_solver=$3
ndom_preprocessor=$4
hostfile=$5
compat=$6
job=$7
srcfile=$8
srcmeth=$9
shift 9 # cannot use $10, $11, ...
restart=$1
postfile=$2
viewfactorsfile=$3
autorst=$4
copy_datfile="-ci $5"
copy_postfile="-cr $6"
scr_dir=$7
dcoup=$8
assem_recov_nthread=$9
shift 9 # cannot use $10, $11, ...
nthread=$1
nsolver=$2
mode=$3
gpu=$4
if [ "$slv" != "" -a "$slv" != "marc" ]; then
slv="-iam sfm"
else
slv=""
fi
if [ "$ndom_fea_solver" != "" -a "$ndom_fea_solver" != "1" ]; then
nprocds="-nprocds $ndom_fea_solver"
else
nprocd=""
if [ "$ndom_preprocessor" != "" -a "$ndom_preprocessor" != "1" ]; then
nprocd="-nprocd $ndom_preprocessor"
else
nprocd=""
fi
fi
if [ "$srcfile" != "" -a "$srcfile" != "-" ]; then
srcfile=`echo $srcfile | sed "s/$SRCEXT$//" | sed "s/$SRCEXTC$//"`
case "$srcmeth" in
-)
srcfile="-u $srcfile"
;;
compsave)
srcfile="-u $srcfile -save y"
;;
runsaved)
srcfile=${srcfile%.*}".marc"
srcfile="-prog $srcfile"
;;
esac
else
srcfile=""
fi
if [ "$restart" != "" -a "$restart" != "-" ]; then
restart=`echo $restart | sed "s/$RSTEXT$//"`
restart="-r $restart"
else
restart=""
fi
if [ "$postfile" != "" -a "$postfile" != "-" ]; then
postfile=`echo $postfile | sed "s/$PSTEXT$//"`
postfile=`echo $postfile | sed "s/$PSTEXTB$//"`
postfile="-pid $postfile"
else
postfile=""
fi
if [ "$viewfactorsfile" != "" -a "$viewfactorsfile" != "-" ]; then
viewfactorsfile=`echo $viewfactorsfile | sed "s/$VWFCEXT$//"`
viewfactorsfile="-vf $viewfactorsfile"
else
viewfactorsfile=""
fi
if [ "$hostfile" != "" -a "$hostfile" != "-" ]; then
hostfile="-ho $hostfile"
else
hostfile=""
fi
if [ "$compat" != "" -a "$compat" != "-" ]; then
compat="-co $compat"
else
compat=""
fi
if [ "$scr_dir" != "" -a "$scr_dir" != "-" ]; then
scr_dir="-sd $scr_dir"
else
scr_dir=""
fi
if [ "$dcoup" != "" -a "$dcoup" != "0" ]; then
dcoup="-dcoup $dcoup"
else
dcoup=""
fi
if [ "$assem_recov_nthread" != "" -a "$assem_recov_nthread" != "1" ]; then
assem_recov_nthread="-nthread_elem $assem_recov_nthread"
else
assem_recov_nthread=""
fi
if [ "$nthread" != "" -a "$nthread" != "0" -a "$nthread" != "1" ]; then
nthread="-nthread $nthread"
else
nthread=""
fi
if [ "$nsolver" != "" -a "$nsolver" != "0" ]; then
nsolver="-nsolver $nsolver"
else
nsolver=""
fi
case "$mode" in
4) mode="-mo i4" ;;
8) mode="-mo i8" ;;
*) mode= ;;
esac
if [ "$gpu" != "" -a "$gpu" != "-" ]; then
gpu="-gpu $gpu"
else
gpu=""
fi
rm -f $job.cnt
rm -f $job.sts
rm -f $job.out
rm -f $job.log
# To prevent a mismatch with the python version used by the solver
# do *not* prepend $MENTAT_INSTALL_DIR/python/bin to environment variable PATH
# unset environment variables PYTHONHOME and PYTHONPATH
unset PYTHONHOME
unset PYTHONPATH
"${DIR}/tools/run_damask" $slv -j $job -v n -b y $nprocds $nprocd -autorst $autorst \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu > /dev/null 2>&1
sleep 1
exit 0

187
installation/mods_MarcMentat/2017/Mentat_bin/submit9
View File

@ -1,187 +0,0 @@
#!/bin/sh
#
# The exit status of this script is read by Mentat.
# Normal exit status is 0.
#
DIR=%INSTALLDIR%/marc%VERSION%
if test $MARCDIR1
then
DIR=$MARCDIR1
fi
if test -z "$DIR"; then
REALCOM="`ls -l $0 |awk '{ print $NF; }'`"
DIRSCRIPT=`dirname $REALCOM`
case $DIRSCRIPT in
\/*)
;;
*)
DIRSCRIPT=`pwd`/$DIRSCRIPT
;;
esac
. $DIRSCRIPT/getarch
DIR="$MENTAT_MARCDIR"
fi
SRCEXT=.f
SRCEXTC=.F
RSTEXT=.t08
PSTEXT=.t19
PSTEXTB=.t16
VWFCEXT=.vfs
slv=$1
version=$2
ndom_fea_solver=$3
ndom_preprocessor=$4
hostfile=$5
compat=$6
job=$7
srcfile=$8
srcmeth=$9
shift 9 # cannot use $10, $11, ...
restart=$1
postfile=$2
viewfactorsfile=$3
autorst=$4
copy_datfile="-ci $5"
copy_postfile="-cr $6"
scr_dir=$7
dcoup=$8
assem_recov_nthread=$9
shift 9 # cannot use $10, $11, ...
nthread=$1
nsolver=$2
mode=$3
gpu=$4
if [ "$slv" != "" -a "$slv" != "marc" ]; then
slv="-iam sfm"
else
slv=""
fi
if [ "$ndom_fea_solver" != "" -a "$ndom_fea_solver" != "1" ]; then
nprocds="-nprocds $ndom_fea_solver"
else
nprocd=""
if [ "$ndom_preprocessor" != "" -a "$ndom_preprocessor" != "1" ]; then
nprocd="-nprocd $ndom_preprocessor"
else
nprocd=""
fi
fi
if [ "$srcfile" != "" -a "$srcfile" != "-" ]; then
srcfile=`echo $srcfile | sed "s/$SRCEXT$//" | sed "s/$SRCEXTC$//"`
case "$srcmeth" in
-)
srcfile="-u $srcfile"
;;
compsave)
srcfile="-u $srcfile -save y"
;;
runsaved)
srcfile=${srcfile%.*}".marc"
srcfile="-prog $srcfile"
;;
esac
else
srcfile=""
fi
if [ "$restart" != "" -a "$restart" != "-" ]; then
restart=`echo $restart | sed "s/$RSTEXT$//"`
restart="-r $restart"
else
restart=""
fi
if [ "$postfile" != "" -a "$postfile" != "-" ]; then
postfile=`echo $postfile | sed "s/$PSTEXT$//"`
postfile=`echo $postfile | sed "s/$PSTEXTB$//"`
postfile="-pid $postfile"
else
postfile=""
fi
if [ "$viewfactorsfile" != "" -a "$viewfactorsfile" != "-" ]; then
viewfactorsfile=`echo $viewfactorsfile | sed "s/$VWFCEXT$//"`
viewfactorsfile="-vf $viewfactorsfile"
else
viewfactorsfile=""
fi
if [ "$hostfile" != "" -a "$hostfile" != "-" ]; then
hostfile="-ho $hostfile"
else
hostfile=""
fi
if [ "$compat" != "" -a "$compat" != "-" ]; then
compat="-co $compat"
else
compat=""
fi
if [ "$scr_dir" != "" -a "$scr_dir" != "-" ]; then
scr_dir="-sd $scr_dir"
else
scr_dir=""
fi
if [ "$dcoup" != "" -a "$dcoup" != "0" ]; then
dcoup="-dcoup $dcoup"
else
dcoup=""
fi
if [ "$assem_recov_nthread" != "" -a "$assem_recov_nthread" != "1" ]; then
assem_recov_nthread="-nthread_elem $assem_recov_nthread"
else
assem_recov_nthread=""
fi
if [ "$nthread" != "" -a "$nthread" != "0" -a "$nthread" != "1" ]; then
nthread="-nthread $nthread"
else
nthread=""
fi
if [ "$nsolver" != "" -a "$nsolver" != "0" ]; then
nsolver="-nsolver $nsolver"
else
nsolver=""
fi
case "$mode" in
4) mode="-mo i4" ;;
8) mode="-mo i8" ;;
*) mode= ;;
esac
if [ "$gpu" != "" -a "$gpu" != "-" ]; then
gpu="-gpu $gpu"
else
gpu=""
fi
rm -f $job.cnt
rm -f $job.sts
rm -f $job.out
rm -f $job.log
# To prevent a mismatch with the python version used by the solver
# do *not* prepend $MENTAT_INSTALL_DIR/python/bin to environment variable PATH
# unset environment variables PYTHONHOME and PYTHONPATH
unset PYTHONHOME
unset PYTHONPATH
"${DIR}/tools/run_damask_l" $slv -j $job -v n -b y $nprocds $nprocd -autorst $autorst \
$srcfile $restart $postfile $viewfactorsfile $hostfile \
$compat $copy_datfile $copy_postfile $scr_dir $dcoup \
$assem_recov_nthread $nthread $nsolver $mode $gpu > /dev/null 2>&1
sleep 1
exit 0

25
installation/mods_MarcMentat/apply_DAMASK_modifications.sh
View File

@ -58,15 +58,9 @@ echo "Editor: $EDITOR"
echo '' echo ''
echo 'adapting Marc tools...' echo 'adapting Marc tools...'
theDIR=$INSTALLDIR/marc$VERSION/tools theDIR=$INSTALLDIR/marc$VERSION/tools
for filename in 'comp_damask' \ for filename in 'comp_damask_mp' \
'comp_damask_l' \
'comp_damask_h' \
'comp_damask_mp' \
'comp_damask_lmp' \ 'comp_damask_lmp' \
'comp_damask_hmp' \ 'comp_damask_hmp' \
'run_damask' \
'run_damask_l' \
'run_damask_h' \
'run_damask_mp' \ 'run_damask_mp' \
'run_damask_lmp' \ 'run_damask_lmp' \
'run_damask_hmp' \ 'run_damask_hmp' \
@ -85,15 +79,9 @@ for filename in 'edit_window' \
'submit4' \ 'submit4' \
'submit5' \ 'submit5' \
'submit6' \ 'submit6' \
'submit7' \
'submit8' \
'submit9' \
'kill4' \ 'kill4' \
'kill5' \ 'kill5' \
'kill6' \ 'kill6'; do
'kill7' \
'kill8' \
'kill9'; do
cp $SCRIPTLOCATION/$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"
@ -122,8 +110,8 @@ echo ''
echo 'setting file access rights...' echo 'setting file access rights...'
for filename in marc$VERSION/tools/run_damask* \ for filename in marc$VERSION/tools/run_damask* \
marc$VERSION/tools/comp_damask* \ marc$VERSION/tools/comp_damask* \
mentat$VERSION/bin/submit{4..9} \ mentat$VERSION/bin/submit{4..6} \
mentat$VERSION/bin/kill{4..9} ; do mentat$VERSION/bin/kill{4..6} ; do
chmod 755 $INSTALLDIR/${filename} chmod 755 $INSTALLDIR/${filename}
done done
@ -142,10 +130,7 @@ if [ -d "$BIN_DIR" ]; then
echo 'creating symlinks ...' echo 'creating symlinks ...'
echo'' echo''
theDIR=$INSTALLDIR/marc$VERSION/tools theDIR=$INSTALLDIR/marc$VERSION/tools
for filename in 'run_damask' \ for filename in 'run_damask_mp' \
'run_damask_l' \
'run_damask_h' \
'run_damask_mp' \
'run_damask_lmp' \ 'run_damask_lmp' \
'run_damask_hmp'; do 'run_damask_hmp'; do
echo ${filename:4}$VERSION echo ${filename:4}$VERSION

12
installation/mods_MarcMentat/installation.txt
View File

@ -21,16 +21,10 @@ The structure of this directory should be (VERSION = 20XX or 20XX.Y)
./installation.txt this text ./installation.txt this text
./apply_MPIE_modifications script file to apply modifications to the installation ./apply_MPIE_modifications script file to apply modifications to the installation
./VERSION/Marc_tools/comp_user.original original file from installation ./VERSION/Marc_tools/comp_user.original original file from installation
./VERSION/Marc_tools/comp_damask modified version using -O1 optimization
./VERSION/Marc_tools/comp_damask_l modified version using -O0 optimization
./VERSION/Marc_tools/comp_damask_h modified version using -O2 optimization
./VERSION/Marc_tools/comp_damask_mp modified version using -O1 optimization and OpenMP ./VERSION/Marc_tools/comp_damask_mp modified version using -O1 optimization and OpenMP
./VERSION/Marc_tools/comp_damask_lmp modified version using -O0 optimization and OpenMP ./VERSION/Marc_tools/comp_damask_lmp modified version using -O0 optimization and OpenMP
./VERSION/Marc_tools/comp_damask_hmp modified version using -O2 optimization and OpenMP ./VERSION/Marc_tools/comp_damask_hmp modified version using -O2 optimization and OpenMP
./VERSION/Marc_tools/run_marc.original original file from installation ./VERSION/Marc_tools/run_marc.original original file from installation
./VERSION/Marc_tools/run_damask modified version using -O1 optimization
./VERSION/Marc_tools/run_damask_l modified version using -O0 optimization
./VERSION/Marc_tools/run_damask_h modified version using -O2 optimization
./VERSION/Marc_tools/run_damask_mp modified version using -O1 optimization and OpenMP ./VERSION/Marc_tools/run_damask_mp modified version using -O1 optimization and OpenMP
./VERSION/Marc_tools/run_damask_lmp modified version using -O0 optimization and OpenMP ./VERSION/Marc_tools/run_damask_lmp modified version using -O0 optimization and OpenMP
./VERSION/Marc_tools/run_damask_hmp modified version using -O2 optimization and OpenMP ./VERSION/Marc_tools/run_damask_hmp modified version using -O2 optimization and OpenMP
@ -42,14 +36,8 @@ The structure of this directory should be (VERSION = 20XX or 20XX.Y)
./VERSION/Mentat_bin/submit4 modified version of original calling run_h_marc ./VERSION/Mentat_bin/submit4 modified version of original calling run_h_marc
./VERSION/Mentat_bin/submit5 modified version of original calling run_marc ./VERSION/Mentat_bin/submit5 modified version of original calling run_marc
./VERSION/Mentat_bin/submit6 modified version of original calling run_l_marc ./VERSION/Mentat_bin/submit6 modified version of original calling run_l_marc
./VERSION/Mentat_bin/submit7 modified version of original calling run_hmp_marc
./VERSION/Mentat_bin/submit8 modified version of original calling run_mp_marc
./VERSION/Mentat_bin/submit9 modified version of original calling run_lmp_marc
./VERSION/Mentat_bin/kill4 kill file for submit4, identical to original kill1 ./VERSION/Mentat_bin/kill4 kill file for submit4, identical to original kill1
./VERSION/Mentat_bin/kill5 kill file for submit5, identical to original kill1 ./VERSION/Mentat_bin/kill5 kill file for submit5, identical to original kill1
./VERSION/Mentat_bin/kill6 kill file for submit6, identical to original kill1 ./VERSION/Mentat_bin/kill6 kill file for submit6, identical to original kill1
./VERSION/Mentat_bin/kill7 kill file for submit7, identical to original kill1
./VERSION/Mentat_bin/kill8 kill file for submit8, identical to original kill1
./VERSION/Mentat_bin/kill9 kill file for submit9, identical to original kill1
./VERSION/Mentat_menus/job_run.ms.original original file from installation ./VERSION/Mentat_menus/job_run.ms.original original file from installation
./VERSION/Mentat_menus/job_run.ms modified version adding DAMASK menu to run menu ./VERSION/Mentat_menus/job_run.ms modified version adding DAMASK menu to run menu

21
src/CMakeLists.txt
View File

@ -17,13 +17,7 @@ list(APPEND OBJECTFILES $<TARGET_OBJECTS:SYSTEM_ROUTINES>)
add_library(PREC OBJECT "prec.f90") add_library(PREC OBJECT "prec.f90")
list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>)
if (PROJECT_NAME STREQUAL "DAMASK_spectral") add_library(DAMASK_INTERFACE OBJECT "DAMASK_interface.f90")
add_library(DAMASK_INTERFACE OBJECT "spectral_interface.f90")
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(DAMASK_INTERFACE OBJECT "FEM_interface.f90")
else ()
message (FATAL_ERROR "Build target (PROJECT_NAME) is not defined")
endif()
add_dependencies(DAMASK_INTERFACE PREC SYSTEM_ROUTINES) add_dependencies(DAMASK_INTERFACE PREC SYSTEM_ROUTINES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>)
@ -57,7 +51,7 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_dependencies(MESH DAMASK_MATH) add_dependencies(MESH DAMASK_MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM") elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEZoo OBJECT "FEZoo.f90") add_library(FEZoo OBJECT "FEM_zoo.f90")
add_dependencies(FEZoo DAMASK_MATH) add_dependencies(FEZoo DAMASK_MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
add_library(MESH OBJECT "meshFEM.f90") add_library(MESH OBJECT "meshFEM.f90")
@ -175,25 +169,24 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral")
"spectral_mech_Basic.f90") "spectral_mech_Basic.f90")
add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES) add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>)
if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES}) add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES})
else() else()
add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90") add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90")
endif() endif()
add_dependencies(DAMASK_spectral SPECTRAL_SOLVER) add_dependencies(DAMASK_spectral SPECTRAL_SOLVER)
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM") elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEM_UTILITIES OBJECT "FEM_utilities.f90") add_library(FEM_UTILITIES OBJECT "FEM_utilities.f90")
add_dependencies(FEM_UTILITIES DAMASK_CPFE) add_dependencies(FEM_UTILITIES DAMASK_CPFE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEM_UTILITIES>)
add_library(FEM_SOLVER OBJECT add_library(FEM_SOLVER OBJECT
"FEM_hydrogenflux.f90"
"FEM_porosity.f90"
"FEM_vacancyflux.f90"
"FEM_damage.f90"
"FEM_thermal.f90"
"FEM_mech.f90") "FEM_mech.f90")
add_dependencies(FEM_SOLVER FEM_UTILITIES) add_dependencies(FEM_SOLVER FEM_UTILITIES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEM_SOLVER>)
add_executable(DAMASK_FEM "DAMASK_FEM_driver.f90") add_executable(DAMASK_FEM "DAMASK_FEM.f90" ${OBJECTFILES})
add_dependencies(DAMASK_FEM FEM_SOLVER) add_dependencies(DAMASK_FEM FEM_SOLVER)
endif() endif()

18
src/CPFEM2.f90
View File

@ -50,8 +50,8 @@ subroutine CPFEM_initAll(el,ip)
IO_init IO_init
use DAMASK_interface use DAMASK_interface
#ifdef FEM #ifdef FEM
use FEZoo, only: & use FEM_Zoo, only: &
FEZoo_init FEM_Zoo_init
#endif #endif
implicit none implicit none
@ -62,7 +62,7 @@ subroutine CPFEM_initAll(el,ip)
call prec_init call prec_init
call IO_init call IO_init
#ifdef FEM #ifdef FEM
call FEZoo_init call FEM_Zoo_init
#endif #endif
call numerics_init call numerics_init
call debug_init call debug_init
@ -196,7 +196,7 @@ end subroutine CPFEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief perform initialization at first call, update variables and call the actual material model !> @brief forwards data after successful increment
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_age() subroutine CPFEM_age()
use prec, only: & use prec, only: &
@ -212,16 +212,6 @@ subroutine CPFEM_age()
debug_levelSelective debug_levelSelective
use FEsolving, only: & use FEsolving, only: &
restartWrite restartWrite
use math, only: &
math_identity2nd, &
math_mul33x33, &
math_det33, &
math_transpose33, &
math_I3, &
math_Mandel3333to66, &
math_Mandel66to3333, &
math_Mandel33to6, &
math_Mandel6to33
use material, only: & use material, only: &
plasticState, & plasticState, &
sourceState, & sourceState, &

654
src/DAMASK_FEM.f90 Normal file
View File

@ -0,0 +1,654 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Driver controlling inner and outer load case looping of the FEM solver
!> @details doing cutbacking, forwarding in case of restart, reporting statistics, writing
!> results
!--------------------------------------------------------------------------------------------------
program DAMASK_FEM
use, intrinsic :: &
iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use prec, only: &
pInt, &
pReal, &
tol_math_check
use DAMASK_interface, only: &
DAMASK_interface_init, &
loadCaseFile, &
getSolverJobName
use IO, only: &
IO_read, &
IO_isBlank, &
IO_open_file, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_error, &
IO_lc, &
IO_intOut, &
IO_warning, &
IO_timeStamp, &
IO_EOF
use debug, only: &
debug_level, &
debug_spectral, &
debug_levelBasic
use math ! need to include the whole module for FFTW
use CPFEM2, only: &
CPFEM_initAll
use FEsolving, only: &
restartWrite, &
restartInc
use numerics, only: &
maxCutBack, &
stagItMax, &
worldrank
use mesh, only: &
mesh_Nboundaries, &
mesh_boundaries, &
geomMesh
use FEM_Utilities, only: &
utilities_init, &
tSolutionState, &
tLoadCase, &
cutBack, &
maxFields, &
nActiveFields, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, &
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID, &
FIELD_MECH_label, &
FIELD_THERMAL_label, &
FIELD_DAMAGE_label, &
FIELD_SOLUTE_label, &
FIELD_MGTWIN_label
use FEM_mech
implicit none
#include <petsc/finclude/petsc.h>
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
integer(pInt), parameter :: FILEUNIT = 234_pInt !< file unit, DAMASK IO does not support newunit feature
integer(pInt), allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing
integer(pInt) :: &
N_def = 0_pInt !< # of rate of deformation specifiers found in load case file
character(len=65536) :: &
line
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
integer(pInt), parameter :: &
subStepFactor = 2_pInt !< for each substep, divide the last time increment by 2.0
real(pReal) :: &
time = 0.0_pReal, & !< elapsed time
time0 = 0.0_pReal, & !< begin of interval
timeinc = 0.0_pReal, & !< current time interval
timeIncOld = 0.0_pReal, & !< previous time interval
remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case
logical :: &
guess !< guess along former trajectory
integer(pInt) :: &
i, &
errorID, &
cutBackLevel = 0_pInt, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0_pInt !< fraction of current time interval
integer(pInt) :: &
currentLoadcase = 0_pInt, & !< current load case
currentFace = 0_pInt, &
inc, & !< current increment in current load case
totalIncsCounter = 0_pInt, & !< total No. of increments
convergedCounter = 0_pInt, & !< No. of converged increments
notConvergedCounter = 0_pInt, & !< No. of non-converged increments
statUnit = 0_pInt, & !< file unit for statistics output
lastRestartWritten = 0_pInt !< total increment No. at which last restart information was written
integer(pInt) :: &
stagIter, &
component
logical :: &
stagIterate
character(len=6) :: loadcase_string
character(len=1024) :: incInfo !< string parsed to solution with information about current load case
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet
PetscInt :: field, dimPlex
PetscErrorCode :: ierr
external :: &
MPI_abort, &
DMGetDimension, &
DMGetLabelSize, &
DMGetLabelIdIS, &
ISDestroy, &
quit
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr)! CHKERRQ(ierr) !< dimension of mesh (2D or 3D)
nActiveFields = 1
allocate(solres(nActiveFields))
!--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases
call IO_open_file(FILEUNIT,trim(loadCaseFile))
rewind(FILEUNIT)
do
line = IO_read(FILEUNIT)
if (trim(line) == IO_EOF) exit
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('$loadcase')
N_def = N_def + 1_pInt
end select
enddo ! count all identifiers to allocate memory and do sanity check
enddo
allocate (loadCases(N_def))
do i = 1, size(loadCases)
allocate(loadCases(i)%fieldBC(nActiveFields))
field = 1
loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID
enddo
do i = 1, size(loadCases)
do field = 1, nActiveFields
select case (loadCases(i)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements
allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents))
end select
do component = 1, loadCases(i)%fieldBC(field)%nComponents
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
enddo
enddo
enddo
!--------------------------------------------------------------------------------------------------
! reading the load case and assign values to the allocated data structure
rewind(FILEUNIT)
do
line = IO_read(FILEUNIT)
if (trim(line) == IO_EOF) exit
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
!--------------------------------------------------------------------------------------------------
! loadcase information
case('$loadcase')
currentLoadCase = IO_intValue(line,chunkPos,i+1_pInt)
case('face')
currentFace = IO_intValue(line,chunkPos,i+1_pInt)
currentFaceSet = -1_pInt
do faceSet = 1, mesh_Nboundaries
if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
enddo
if (currentFaceSet < 0_pInt) call IO_error(error_ID = errorID, ext_msg = 'invalid BC')
case('t','time','delta') ! increment time
loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1_pInt)
case('n','incs','increments','steps') ! number of increments
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt)
loadCases(currentLoadCase)%logscale = 1_pInt
case('freq','frequency','outputfreq') ! frequency of result writings
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt)
case('r','restart','restartwrite') ! frequency of writing restart information
loadCases(currentLoadCase)%restartfrequency = &
max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt))
case('guessreset','dropguessing')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
!--------------------------------------------------------------------------------------------------
! boundary condition information
case('x') ! X displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_X_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('y') ! Y displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Y_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('z') ! Z displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Z_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('temp','temperature') ! thermal field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_THERMAL_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_THERMAL_T_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('mgtwin') ! mgtwin field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MGTWIN_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MGTWIN_PHI_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('damage')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_DAMAGE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_DAMAGE_PHI_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cv')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CV_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('ch')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CH_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('chpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CHPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvah')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaH_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvahpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaHPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
end select
enddo; enddo
close(FILEUNIT)
!--------------------------------------------------------------------------------------------------
! consistency checks and output of load case
loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase
errorID = 0_pInt
if (worldrank == 0) then
checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases)
write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
write(6,'(2x,a)') 'drop guessing along trajectory'
do field = 1_pInt, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label)
case(FIELD_THERMAL_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_THERMAL_label)
case(FIELD_DAMAGE_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_DAMAGE_label)
case(FIELD_MGTWIN_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MGTWIN_label)
case(FIELD_SOLUTE_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_SOLUTE_label)
end select
do faceSet = 1_pInt, mesh_Nboundaries
do component = 1_pInt, loadCases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) &
write(6,'(4x,a,i2,a,i2,a,f12.7)') 'Face ', mesh_boundaries(faceSet), &
' Component ', component, &
' Value ', loadCases(currentLoadCase)%fieldBC(field)% &
componentBC(component)%Value(faceSet)
enddo
enddo
enddo
write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count
write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs
if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', &
loadCases(currentLoadCase)%outputfrequency
write(6,'(2x,a,i5,/)') 'restart frequency: ', &
loadCases(currentLoadCase)%restartfrequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases
endif
!--------------------------------------------------------------------------------------------------
! doing initialization depending on selected solver
call Utilities_init()
do field = 1, nActiveFields
select case (loadCases(1)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_init(loadCases(1)%fieldBC(field))
end select
enddo
!--------------------------------------------------------------------------------------------------
! loopping over loadcases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
time0 = time ! currentLoadCase start time
if (loadCases(currentLoadCase)%followFormerTrajectory) then
guess = .true.
else
guess = .false. ! change of load case, homogeneous guess for the first inc
endif
!--------------------------------------------------------------------------------------------------
! loop oper incs defined in input file for current currentLoadCase
incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt
!--------------------------------------------------------------------------------------------------
! forwarding time
timeIncOld = timeinc
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/loadCases(currentLoadCase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
else
if (currentLoadCase == 1_pInt) then ! 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
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))
endif
else ! not-1st currentLoadCase of logarithmic scale
timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal)))
endif
endif
timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step
forwarding: if(totalIncsCounter >= restartInc) then
stepFraction = 0_pInt
!--------------------------------------------------------------------------------------------------
! loop over sub incs
subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt)
time = time + timeinc ! forward time
stepFraction = stepFraction + 1_pInt
remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc
!--------------------------------------------------------------------------------------------------
! report begin of new increment
if (worldrank == 0) then
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases)
flush(6)
write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel
endif
!--------------------------------------------------------------------------------------------------
! forward fields
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
enddo
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0_pInt
stagIterate = .true.
do while (stagIterate)
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
solres(field) = FEM_mech_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
if(.not. solres(field)%converged) exit ! no solution found
enddo
stagIter = stagIter + 1_pInt
stagIterate = stagIter < stagItMax .and. &
all(solres(:)%converged) .and. &
.not. all(solres(:)%stagConverged)
enddo
! check solution
cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back
if (worldrank == 0) &
write(6,'(/,a)') ' cut back detected'
cutBack = .True.
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt
time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pReal
else ! default behavior, exit if spectral solver does not converge
call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written (e.g. for regridding) ! continue from non-converged solution and start guessing after accepted (sub)inc
endif
else
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
endif
if (.not. cutBack) then
if (worldrank == 0) write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
endif
enddo subIncLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
if(all(solres(:)%converged)) then ! report converged inc
convergedCounter = convergedCounter + 1_pInt
if (worldrank == 0) then
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &
' increment ', totalIncsCounter, ' converged'
endif
else
if (worldrank == 0) then
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif
notConvergedCounter = notConvergedCounter + 1_pInt
endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
if (worldrank == 0) then
write(6,'(1/,a)') ' ... writing results to file ......................................'
endif
endif
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 ! ToDo first call to CPFEM_general will write?
restartWrite = .true.
lastRestartWritten = inc
endif
else forwarding
time = time + timeinc
guess = .true.
endif forwarding
enddo incLooping
enddo loadCaseLooping
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
if (worldrank == 0) then
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
' %) increments converged!'
endif
if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;)
end program DAMASK_FEM
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief quit subroutine to mimic behavior of FEM solvers
!> @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
!> 2 signals request for regridding, increment of last saved restart information is written to
!> stderr. Exit code 3 signals no severe problems, but some increments did not converge
!--------------------------------------------------------------------------------------------------
subroutine quit(stop_id)
use prec, only: &
pInt
implicit none
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id < 0_pInt) then ! trigger regridding
write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
stop 2
endif
if (stop_id == 3_pInt) stop 3 ! not all incs converged
stop 1 ! error (message from IO_error)
end subroutine quit

160
src/spectral_interface.f90 → src/DAMASK_interface.f90
View File

@ -1,9 +1,11 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Jaeyong Jung, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Interfacing between the spectral solver and the material subroutines provided !> @brief Interfacing between the PETSc-based solvers and the material subroutines provided
!! by DAMASK !! by DAMASK
!> @details Interfacing between the spectral solver and the material subroutines provided !> @details Interfacing between the PETSc-based solvers and the material subroutines provided
!> by DAMASK. Interpretating the command line arguments to get load case, geometry file, !> by DAMASK. Interpretating the command line arguments to get load case, geometry file,
!> and working directory. !> and working directory.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -13,12 +15,11 @@ module DAMASK_interface
implicit none implicit none
private private
logical, public, protected :: appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding) integer(pInt), public, protected :: &
integer(pInt), public, protected :: spectralRestartInc = 0_pInt !< Increment at which calculation starts interface_restartInc = 0_pInt !< Increment at which calculation starts
character(len=1024), public, protected :: & character(len=1024), public, protected :: &
geometryFile = '', & !< parameter given for geometry file geometryFile = '', & !< parameter given for geometry file
loadCaseFile = '' !< parameter given for load case file loadCaseFile = '' !< parameter given for load case file
character(len=1024), private :: workingDirectory
public :: & public :: &
getSolverJobName, & getSolverJobName, &
@ -44,22 +45,36 @@ subroutine DAMASK_interface_init()
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
#if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR!=9 #if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR!=9
=================================================================================================== ===================================================================================================
3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x
===================================================================================================
======= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ===========================================
========== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ========================================
============= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =====================================
================ THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ==================================
=================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ===============================
====================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ============================
========================= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ========================= ========================= THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =========================
============================ THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ======================
=============================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ===================
================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ================
===================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x =============
======================================== THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x ==========
===================================================================================================
3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x
=================================================================================================== ===================================================================================================
#endif #endif
use PETScSys use PETScSys
use system_routines, only: & use system_routines, only: &
getHostName getHostName, &
getCWD
implicit none implicit none
character(len=1024) :: & character(len=1024) :: &
commandLine, & !< command line call as string commandLine, & !< command line call as string
loadcaseArg = '', & !< -l argument given to DAMASK_spectral.exe loadcaseArg = '', & !< -l argument given to the executable
geometryArg = '', & !< -g argument given to DAMASK_spectral.exe geometryArg = '', & !< -g argument given to the executable
workingDirArg = '', & !< -w argument given to DAMASK_spectral.exe workingDirArg = '', & !< -w argument given to the executable
hostName, & !< name of machine on which DAMASK_spectral.exe is execute (might require export HOSTNAME) userName !< name of user calling the executable
userName, & !< name of user calling DAMASK_spectral.exe
tag
integer :: & integer :: &
i, & i, &
#ifdef _OPENMP #ifdef _OPENMP
@ -72,7 +87,6 @@ subroutine DAMASK_interface_init()
integer, dimension(8) :: & integer, dimension(8) :: &
dateAndTime ! type default integer dateAndTime ! type default integer
PetscErrorCode :: ierr PetscErrorCode :: ierr
logical :: error
external :: & external :: &
quit,& quit,&
PETScErrorF, & ! is called in the CHKERRQ macro PETScErrorF, & ! is called in the CHKERRQ macro
@ -110,7 +124,7 @@ subroutine DAMASK_interface_init()
endif mainProcess endif mainProcess
call date_and_time(values = dateAndTime) call date_and_time(values = dateAndTime)
write(6,'(/,a)') ' <<<+- DAMASK_spectral -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
write(6,'(a,/)') ' Roters et al., Computational Materials Science, 2018' write(6,'(a,/)') ' Roters et al., Computational Materials Science, 2018'
write(6,'(/,a)') ' Version: '//DAMASKVERSION write(6,'(/,a)') ' Version: '//DAMASKVERSION
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',& write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',&
@ -120,7 +134,6 @@ subroutine DAMASK_interface_init()
dateAndTime(6),':',& dateAndTime(6),':',&
dateAndTime(7) dateAndTime(7)
write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
#include "compilation_info.f90" #include "compilation_info.f90"
call get_command(commandLine) call get_command(commandLine)
@ -129,9 +142,8 @@ subroutine DAMASK_interface_init()
select case(IIO_stringValue(commandLine,chunkPos,i)) ! extract key select case(IIO_stringValue(commandLine,chunkPos,i)) ! extract key
case ('-h','--help') case ('-h','--help')
write(6,'(a)') ' #######################################################################' write(6,'(a)') ' #######################################################################'
write(6,'(a)') ' DAMASK_spectral:' write(6,'(a)') ' DAMASK Command Line Interface:'
write(6,'(a)') ' The spectral method boundary value problem solver for' write(6,'(a)') ' For PETSc-based solvers for the Düsseldorf Advanced Material Simulation Kit'
write(6,'(a)') ' the Düsseldorf Advanced Material Simulation Kit'
write(6,'(a,/)')' #######################################################################' write(6,'(a,/)')' #######################################################################'
write(6,'(a,/)')' Valid command line switches:' write(6,'(a,/)')' Valid command line switches:'
write(6,'(a)') ' --geom (-g, --geometry)' write(6,'(a)') ' --geom (-g, --geometry)'
@ -141,23 +153,14 @@ subroutine DAMASK_interface_init()
write(6,'(a)') ' --help (-h)' write(6,'(a)') ' --help (-h)'
write(6,'(/,a)')' -----------------------------------------------------------------------' write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Mandatory arguments:' write(6,'(a)') ' Mandatory arguments:'
write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom.geom' write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom'
write(6,'(a)') ' Specifies the location of the geometry definition file,' write(6,'(a)') ' Specifies the location of the geometry definition file.'
write(6,'(a)') ' if no extension is given, .geom will be appended.' write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile'
write(6,'(a)') ' "PathToGeomFile" will be the working directory if not specified' write(6,'(a)') ' Specifies the location of the load case definition file.'
write(6,'(a)') ' via --workingdir.'
write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' directory.'
write(6,'(a)') ' For further configuration place "numerics.config"'
write(6,'(a)')' and "numerics.config" in that directory.'
write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile.load'
write(6,'(a)') ' Specifies the location of the load case definition file,'
write(6,'(a)') ' if no extension is given, .load will be appended.'
write(6,'(/,a)')' -----------------------------------------------------------------------' write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Optional arguments:' write(6,'(a)') ' Optional arguments:'
write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory' write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory'
write(6,'(a)') ' Specifies the working directory and overwrites the default' write(6,'(a)') ' Specifies the working directory and overwrites the default ./'
write(6,'(a)') ' "PathToGeomFile".'
write(6,'(a)') ' Make sure the file "material.config" exists in the working' write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' directory.' write(6,'(a)') ' directory.'
write(6,'(a)') ' For further configuration place "numerics.config"' write(6,'(a)') ' For further configuration place "numerics.config"'
@ -166,7 +169,7 @@ subroutine DAMASK_interface_init()
write(6,'(a)') ' Reads in increment XX and continues with calculating' write(6,'(a)') ' Reads in increment XX and continues with calculating'
write(6,'(a)') ' increment XX+1 based on this.' write(6,'(a)') ' increment XX+1 based on this.'
write(6,'(a)') ' Appends to existing results file' write(6,'(a)') ' Appends to existing results file'
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' write(6,'(a)') ' "NameOfGeom_NameOfLoadFile".'
write(6,'(a)') ' Works only if the restart information for increment XX' write(6,'(a)') ' Works only if the restart information for increment XX'
write(6,'(a)') ' is available in the working directory.' write(6,'(a)') ' is available in the working directory.'
write(6,'(/,a)')' -----------------------------------------------------------------------' write(6,'(/,a)')' -----------------------------------------------------------------------'
@ -182,8 +185,7 @@ subroutine DAMASK_interface_init()
if (i < chunkPos(1)) workingDirArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt)) if (i < chunkPos(1)) workingDirArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt))
case ('-r', '--rs', '--restart') case ('-r', '--rs', '--restart')
if (i < chunkPos(1)) then if (i < chunkPos(1)) then
spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt) interface_restartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
appendToOutFile = .true.
endif endif
end select end select
enddo enddo
@ -193,26 +195,25 @@ subroutine DAMASK_interface_init()
call quit(1_pInt) call quit(1_pInt)
endif endif
workingDirectory = trim(setWorkingDirectory(trim(workingDirArg))) if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg))
geometryFile = getGeometryFile(geometryArg) geometryFile = getGeometryFile(geometryArg)
loadCaseFile = getLoadCaseFile(loadCaseArg) loadCaseFile = getLoadCaseFile(loadCaseArg)
call get_environment_variable('USER',userName) call get_environment_variable('USER',userName)
error = getHostName(hostName) ! ToDo: https://stackoverflow.com/questions/8953424/how-to-get-the-username-in-c-c-in-linux
write(6,'(a,a)') ' Host name: ', trim(hostName) write(6,'(a,a)') ' Host name: ', trim(getHostName())
write(6,'(a,a)') ' User name: ', trim(userName) write(6,'(a,a)') ' User name: ', trim(userName)
write(6,'(a,a)') ' Command line call: ', trim(commandLine) write(6,'(a,a)') ' Command line call: ', trim(commandLine)
if (len(trim(workingDirArg)) > 0) & if (len(trim(workingDirArg)) > 0) &
write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg)
write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg)
write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg) write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg)
write(6,'(a,a)') ' Working directory: ', trim(workingDirectory) write(6,'(a,a)') ' Working directory: ', trim(getCWD())
write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) write(6,'(a,a)') ' Geometry file: ', trim(geometryFile)
write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile)
write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName())
if (SpectralRestartInc > 0_pInt) & if (interface_restartInc > 0_pInt) &
write(6,'(a,i6.6)') ' Restart from increment: ', spectralRestartInc write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc
write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile
end subroutine DAMASK_interface_init end subroutine DAMASK_interface_init
@ -221,38 +222,32 @@ end subroutine DAMASK_interface_init
!> @brief extract working directory from given argument or from location of geometry file, !> @brief extract working directory from given argument or from location of geometry file,
!! possibly converting relative arguments to absolut path !! possibly converting relative arguments to absolut path
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
character(len=1024) function setWorkingDirectory(workingDirectoryArg) subroutine setWorkingDirectory(workingDirectoryArg)
use system_routines, only: & use system_routines, only: &
getCWD, & getCWD, &
setCWD setCWD
implicit none implicit none
character(len=*), intent(in) :: workingDirectoryArg !< working directory argument character(len=*), intent(in) :: workingDirectoryArg !< working directory argument
logical :: error character(len=1024) :: workingDirectory !< working directory argument
external :: quit external :: quit
logical :: error
wdGiven: if (len(workingDirectoryArg)>0) then absolutePath: if (workingDirectoryArg(1:1) == '/') then
absolutePath: if (workingDirectoryArg(1:1) == '/') then workingDirectory = workingDirectoryArg
setWorkingDirectory = workingDirectoryArg else absolutePath
else absolutePath workingDirectory = getCWD()
error = getCWD(setWorkingDirectory) workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg
if (error) call quit(1_pInt) endif absolutePath
setWorkingDirectory = trim(setWorkingDirectory)//'/'//workingDirectoryArg
endif absolutePath
else wdGiven
error = getCWD(setWorkingDirectory) ! relative path given as command line argument
if (error) call quit(1_pInt)
endif wdGiven
setWorkingDirectory = trim(rectifyPath(setWorkingDirectory)) workingDirectory = trim(rectifyPath(workingDirectory))
error = setCWD(trim(workingDirectory))
error = setCWD(trim(setWorkingDirectory))
if(error) then if(error) then
write(6,'(a20,a,a16)') ' working directory "',trim(setWorkingDirectory),'" does not exist' write(6,'(a20,a,a16)') ' working directory "',trim(workingDirectory),'" does not exist'
call quit(1_pInt) call quit(1_pInt)
endif endif
end function setWorkingDirectory end subroutine setWorkingDirectory
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -284,22 +279,15 @@ end function getSolverJobName
!> @brief basename of geometry file with extension from command line arguments !> @brief basename of geometry file with extension from command line arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
character(len=1024) function getGeometryFile(geometryParameter) character(len=1024) function getGeometryFile(geometryParameter)
use system_routines, only: &
getCWD
implicit none implicit none
character(len=1024), intent(in) :: & character(len=1024), intent(in) :: geometryParameter
geometryParameter
integer :: posExt, posSep
external :: quit
getGeometryFile = trim(geometryParameter) getGeometryFile = trim(geometryParameter)
posExt = scan(getGeometryFile,'.',back=.true.) if (scan(getGeometryFile,'/') /= 1) getGeometryFile = trim(getCWD())//'/'//trim(getGeometryFile)
posSep = scan(getGeometryFile,'/',back=.true.) getGeometryFile = makeRelativePath(trim(getCWD()), getGeometryFile)
if (posExt <= posSep) getGeometryFile = trim(getGeometryFile)//('.geom')
if (scan(getGeometryFile,'/') /= 1) &
getGeometryFile = trim(workingDirectory)//'/'//trim(getGeometryFile)
getGeometryFile = makeRelativePath(workingDirectory, getGeometryFile)
end function getGeometryFile end function getGeometryFile
@ -309,22 +297,15 @@ end function getGeometryFile
!> @brief relative path of loadcase from command line arguments !> @brief relative path of loadcase from command line arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
character(len=1024) function getLoadCaseFile(loadCaseParameter) character(len=1024) function getLoadCaseFile(loadCaseParameter)
use system_routines, only: &
getCWD
implicit none implicit none
character(len=1024), intent(in) :: & character(len=1024), intent(in) :: loadCaseParameter
loadCaseParameter
integer :: posExt, posSep
external :: quit
getLoadCaseFile = trim(loadCaseParameter) getLoadCaseFile = trim(loadCaseParameter)
posExt = scan(getLoadCaseFile,'.',back=.true.) if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = trim(getCWD())//'/'//trim(getLoadCaseFile)
posSep = scan(getLoadCaseFile,'/',back=.true.) getLoadCaseFile = makeRelativePath(trim(getCWD()), getLoadCaseFile)
if (posExt <= posSep) getLoadCaseFile = trim(getLoadCaseFile)//('.load')
if (scan(getLoadCaseFile,'/') /= 1) &
getLoadCaseFile = trim(workingDirectory)//'/'//trim(getLoadCaseFile)
getLoadCaseFile = makeRelativePath(workingDirectory, getLoadCaseFile)
end function getLoadCaseFile end function getLoadCaseFile
@ -337,21 +318,20 @@ function rectifyPath(path)
implicit none implicit none
character(len=*) :: path character(len=*) :: path
character(len=len_trim(path)) :: rectifyPath character(len=1024) :: rectifyPath
integer :: i,j,k,l ! no pInt integer :: i,j,k,l ! no pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! remove /./ from path ! remove /./ from path
l = len_trim(path) rectifyPath = trim(path)
rectifyPath = path l = len_trim(rectifyPath)
do i = l,3,-1 do i = l,3,-1
if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! remove // from path ! remove // from path
l = len_trim(path) l = len_trim(rectifyPath)
rectifyPath = path
do i = l,2,-1 do i = l,2,-1
if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' ' if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' '
enddo enddo

15
src/DAMASK_spectral.f90
View File

@ -20,14 +20,12 @@ program DAMASK_spectral
pReal, & pReal, &
tol_math_check, & tol_math_check, &
dNeq dNeq
use system_routines, only: &
getCWD
use DAMASK_interface, only: & use DAMASK_interface, only: &
DAMASK_interface_init, & DAMASK_interface_init, &
loadCaseFile, & loadCaseFile, &
geometryFile, & geometryFile, &
getSolverJobName, & getSolverJobName, &
appendToOutFile interface_restartInc
use IO, only: & use IO, only: &
IO_read, & IO_read, &
IO_isBlank, & IO_isBlank, &
@ -386,8 +384,7 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write header of output file ! write header of output file
if (worldrank == 0) then if (worldrank == 0) then
if (.not. appendToOutFile) then ! after restart, append to existing results file writeHeader: if (interface_restartInc < 1_pInt) then
if (getCWD(workingDir)) call IO_error(106_pInt,ext_msg=trim(workingDir))
open(newunit=resUnit,file=trim(getSolverJobName())//& open(newunit=resUnit,file=trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED',status='REPLACE') '.spectralOut',form='UNFORMATTED',status='REPLACE')
write(resUnit) 'load:', trim(loadCaseFile) ! ... and write header write(resUnit) 'load:', trim(loadCaseFile) ! ... and write header
@ -410,10 +407,10 @@ program DAMASK_spectral
if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) & if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) &
write(6,'(/,a)') ' header of result and statistics file written out' write(6,'(/,a)') ' header of result and statistics file written out'
flush(6) flush(6)
else ! open new files ... else writeHeader
open(newunit=statUnit,file=trim(getSolverJobName())//& open(newunit=statUnit,file=trim(getSolverJobName())//&
'.sta',form='FORMATTED', position='APPEND', status='OLD') '.sta',form='FORMATTED', position='APPEND', status='OLD')
endif endif writeHeader
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -434,7 +431,7 @@ program DAMASK_spectral
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(error_ID=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 writeUndeformed: if (interface_restartInc < 1_pInt) then
write(6,'(1/,a)') ' ... writing initial configuration to file ........................' write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
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, & ! QUESTION: why not starting i at 0 instead of murky 1? outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & ! QUESTION: why not starting i at 0 instead of murky 1?
@ -446,7 +443,7 @@ program DAMASK_spectral
if (ierr /= 0_pInt) call IO_error(error_ID=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
endif endif writeUndeformed
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! looping over loadcases ! looping over loadcases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)

734
src/FEM_mech.f90 Normal file
View File

@ -0,0 +1,734 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief FEM PETSc solver
!--------------------------------------------------------------------------------------------------
module FEM_mech
#include <petsc/finclude/petsc.h>
use PETScdmda
use PETScsnes
use PETScDM
use PETScDMplex
use prec, only: &
pInt, &
pReal
use math, only: &
math_I3
use FEM_utilities, only: &
tSolutionState, &
tFieldBC, &
tComponentBC
use numerics, only: &
worldrank, &
worldsize
use mesh, only: &
mesh_Nboundaries, &
mesh_boundaries
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type tSolutionParams
type(tFieldBC) :: fieldBC
real(pReal) :: timeinc
real(pReal) :: timeincOld
end type tSolutionParams
type(tSolutionParams), private :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
SNES, private :: mech_snes
Vec, private :: solution, solution_rate, solution_local
PetscInt, private :: dimPlex, cellDof, nQuadrature, nBasis
PetscReal, allocatable, target,dimension(:), private :: qPoints, qWeights
MatNullSpace, private :: matnull
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
character(len=1024), private :: incInfo
real(pReal), private, dimension(3,3) :: &
P_av = 0.0_pReal
logical, private :: ForwardData
real(pReal), parameter, private :: eps = 1.0e-18_pReal
public :: &
FEM_mech_init, &
FEM_mech_solution ,&
FEM_mech_forward, &
FEM_mech_destroy
external :: &
MatZeroRowsColumnsLocalIS, &
PetscQuadratureCreate, &
PetscFECreateDefault, &
PetscFESetQuadrature, &
PetscFEGetDimension, &
PetscFEDestroy, &
PetscFEGetDualSpace, &
PetscQuadratureDestroy, &
PetscDSSetDiscretization, &
PetscDSGetTotalDimension, &
PetscDSGetDiscretization, &
PetscDualSpaceGetFunctional, &
DMGetLabelSize, &
DMSNESSetFunctionLocal, &
DMSNESSetJacobianLocal, &
SNESSetOptionsPrefix, &
SNESSetConvergenceTest, &
PetscObjectSetName
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_init(fieldBC)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO, only: &
IO_timeStamp, &
IO_error
use DAMASK_interface, only: &
getSolverJobName
use mesh, only: &
geomMesh
use numerics, only: &
worldrank, &
itmax, &
integrationOrder
use FEM_Zoo, only: &
FEM_Zoo_nQuadrature, &
FEM_Zoo_QuadraturePoints, &
FEM_Zoo_QuadratureWeights
implicit none
type(tFieldBC), intent(in) :: fieldBC
DM :: mech_mesh
PetscFE :: mechFE
PetscQuadrature :: mechQuad, functional
PetscDS :: mechDS
PetscDualSpace :: mechDualSpace
DMLabel :: BCLabel
PetscInt, allocatable, target :: numComp(:), numDoF(:), bcField(:)
PetscInt, pointer :: pNumComp(:), pNumDof(:), pBcField(:), pBcPoint(:)
PetscInt :: numBC, bcSize
IS :: bcPoint
IS, allocatable, target :: bcComps(:), bcPoints(:)
IS, pointer :: pBcComps(:), pBcPoints(:)
PetscSection :: section
PetscInt :: field, faceSet, topologDim, nNodalPoints
PetscReal, pointer :: qPointsP(:), qWeightsP(:), &
nodalPointsP(:), nodalWeightsP(:)
PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:)
PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:)
PetscReal :: detJ
PetscReal, allocatable, target :: v0(:), cellJ(:), invcellJ(:), cellJMat(:,:)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscInt :: cellStart, cellEnd, cell, basis
character(len=7) :: prefix = 'mechFE_'
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! Setup FEM mech mesh
call DMClone(geomMesh,mech_mesh,ierr); CHKERRQ(ierr)
call DMGetDimension(mech_mesh,dimPlex,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! Setup FEM mech discretization
allocate(qPoints(dimPlex*FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
allocate(qWeights(FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
qPoints = FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p
qWeights = FEM_Zoo_QuadratureWeights(dimPlex,integrationOrder)%p
nQuadrature = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
qPointsP => qPoints
qWeightsP => qWeights
call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr)
call PetscQuadratureSetData(mechQuad,dimPlex,nQuadrature,qPointsP,qWeightsP,ierr)
CHKERRQ(ierr)
call PetscFECreateDefault(mech_mesh,dimPlex,dimPlex,PETSC_TRUE,prefix, &
integrationOrder,mechFE,ierr); CHKERRQ(ierr)
call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr)
call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr)
call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr)
call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! Setup FEM mech boundary conditions
call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr)
call DMGetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr)
allocate(numComp(1), source=dimPlex); pNumComp => numComp
allocate(numDof(dimPlex+1), source = 0); pNumDof => numDof
do topologDim = 0, dimPlex
call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr)
CHKERRQ(ierr)
call PetscSectionGetDof(section,cellStart,numDof(topologDim+1),ierr)
CHKERRQ(ierr)
enddo
numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1
enddo; enddo
allocate(bcField(numBC), source=0); pBcField => bcField
allocate(bcComps(numBC)); pBcComps => bcComps
allocate(bcPoints(numBC)); pBcPoints => bcPoints
numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) then
numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,bcComps(numBC),ierr)
CHKERRQ(ierr)
call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr)
CHKERRQ(ierr)
if (bcSize > 0) then
call DMGetStratumIS(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr)
CHKERRQ(ierr)
call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,bcPoints(numBC),ierr)
CHKERRQ(ierr)
call ISRestoreIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISDestroy(bcPoint,ierr); CHKERRQ(ierr)
else
call ISCreateGeneral(PETSC_COMM_WORLD,0,[0],PETSC_COPY_VALUES,bcPoints(numBC),ierr)
CHKERRQ(ierr)
endif
endif
enddo; enddo
call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, &
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS, &
section,ierr)
CHKERRQ(ierr)
call DMSetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr)
do faceSet = 1, numBC
call ISDestroy(bcPoints(faceSet),ierr); CHKERRQ(ierr)
enddo
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr);CHKERRQ(ierr)
call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr)
call SNESSetDM(mech_snes,mech_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver
call DMCreateGlobalVector(mech_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs
call DMCreateGlobalVector(mech_mesh,solution_rate ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMCreateLocalVector (mech_mesh,solution_local ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMSNESSetFunctionLocal(mech_mesh,FEM_mech_formResidual,PETSC_NULL_VEC,ierr) !< function to evaluate residual forces
CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mech_mesh,FEM_mech_formJacobian,PETSC_NULL_VEC,ierr) !< function to evaluate stiffness matrix
CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures
call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr)
call SNESSetTolerances(mech_snes,1.0,0.0,0.0,itmax,itmax,ierr)
CHKERRQ(ierr)
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! init fields
call VecSet(solution ,0.0,ierr); CHKERRQ(ierr)
call VecSet(solution_rate ,0.0,ierr); CHKERRQ(ierr)
allocate(x_scal(cellDof))
allocate(nodalPoints (dimPlex))
allocate(nodalWeights(1))
nodalPointsP => nodalPoints
nodalWeightsP => nodalWeights
allocate(v0(dimPlex))
allocate(cellJ(dimPlex*dimPlex))
allocate(invcellJ(dimPlex*dimPlex))
allocate(cellJMat(dimPlex,dimPlex))
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMGetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSGetDiscretization(mechDS,0,mechFE,ierr)
CHKERRQ(ierr)
call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr)
call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
do basis = 0, nBasis-1
call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr)
CHKERRQ(ierr)
call PetscQuadratureGetData(functional,dimPlex,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
CHKERRQ(ierr)
x_scal(basis*dimPlex+1:(basis+1)*dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0)
enddo
px_scal => x_scal
call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr)
CHKERRQ(ierr)
enddo
end subroutine FEM_mech_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM load step
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function FEM_mech_solution( &
incInfoIn,timeinc,timeinc_old,fieldBC)
use numerics, only: &
itmax
use FEsolving, only: &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment
type(tFieldBC), intent(in) :: &
fieldBC
character(len=*), intent(in) :: &
incInfoIn
!--------------------------------------------------------------------------------------------------
!
PetscErrorCode :: ierr
SNESConvergedReason :: reason
incInfo = incInfoIn
FEM_mech_solution%converged =.false.
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%timeinc = timeinc
params%timeincOld = timeinc_old
params%fieldBC = fieldBC
call SNESSolve(mech_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mech_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
FEM_mech_solution%converged = .false.
FEM_mech_solution%iterationsNeeded = itmax
else ! >= 1 proper convergence (or terminally ill)
FEM_mech_solution%converged = .true.
call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr)
CHKERRQ(ierr)
endif
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end function FEM_mech_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the FEM residual vector
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
use numerics, only: &
BBarStabilisation
use FEM_utilities, only: &
utilities_projectBCValues, &
utilities_constitutiveResponse
use homogenization, only: &
materialpoint_F, &
materialpoint_P
use math, only: &
math_det33, &
math_inv33
use FEsolving, only: &
terminallyIll
implicit none
DM :: dm_local
PetscDS :: prob
Vec :: x_local, f_local, xx_local
PetscSection :: section
PetscScalar, dimension(:), pointer :: x_scal, pf_scal
PetscScalar, target :: f_scal(cellDof)
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), &
invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal, pointer :: basisField(:), basisFieldDer(:)
PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx
PetscReal :: detFAvg
PetscReal :: BMat(dimPlex*dimPlex,cellDof)
PetscObject :: dummy
PetscInt :: bcSize
IS :: bcPoints
PetscErrorCode :: ierr
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
CHKERRQ(ierr)
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
enddo; enddo
!--------------------------------------------------------------------------------------------------
! evaluate field derivatives
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
do qPt = 0, nQuadrature-1
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx )*dimPlex+1: &
(qPt*nBasis*dimPlex+cidx+1)*dimPlex ))
enddo
enddo
materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo
if (BBarStabilisation) then
detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature))
do qPt = 1, nQuadrature
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = &
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1)* &
(detFAvg/math_det33(materialpoint_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex))
enddo
endif
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr)
enddo
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
ForwardData = .false.
!--------------------------------------------------------------------------------------------------
! integrating residual
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
f_scal = 0.0
do qPt = 0, nQuadrature-1
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx )*dimPlex+1: &
(qPt*nBasis*dimPlex+cidx+1)*dimPlex ))
enddo
enddo
f_scal = f_scal + &
matmul(transpose(BMat), &
reshape(transpose(materialpoint_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), &
shape=[dimPlex*dimPlex]))*qWeights(qPt+1)
enddo
f_scal = f_scal*abs(detJ)
pf_scal => f_scal
call DMPlexVecSetClosure(dm_local,section,f_local,cell,pf_scal,ADD_VALUES,ierr)
CHKERRQ(ierr)
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr)
enddo
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_formResidual
!--------------------------------------------------------------------------------------------------
!> @brief forms the FEM stiffness matrix
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
use numerics, only: &
BBarStabilisation
use homogenization, only: &
materialpoint_dPdF, &
materialpoint_F
use math, only: &
math_inv33, &
math_identity2nd, &
math_det33
use FEM_utilities, only: &
utilities_projectBCValues
implicit none
DM :: dm_local
PetscDS :: prob
Vec :: x_local, xx_local
Mat :: Jac_pre, Jac
PetscSection :: section, gSection
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), &
invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal, dimension(:), pointer :: basisField, basisFieldDer
PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx
PetscScalar, target :: K_e (cellDof,cellDof), &
K_eA (cellDof,cellDof), &
K_eB (cellDof,cellDof), &
K_eVec(cellDof*cellDof)
PetscReal :: BMat (dimPlex*dimPlex,cellDof), &
BMatAvg(dimPlex*dimPlex,cellDof), &
MatA (dimPlex*dimPlex,cellDof), &
MatB (1 ,cellDof)
PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscReal, dimension(3,3) :: F = math_I3, FAvg, FInv
PetscObject :: dummy
PetscInt :: bcSize
IS :: bcPoints
PetscErrorCode :: ierr
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr)
call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr)
call MatZeroEntries(Jac,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetDefaultGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
enddo; enddo
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ, shape = [dimPlex,dimPlex])
K_eA = 0.0
K_eB = 0.0
MatB = 0.0
FAvg = 0.0
BMatAvg = 0.0
do qPt = 0, nQuadrature-1
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx )*dimPlex+1: &
(qPt*nBasis*dimPlex+cidx+1)*dimPlex ))
enddo
enddo
MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
if (BBarStabilisation) then
F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex])
FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
K_eB = K_eB - &
matmul(transpose(matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex*dimPlex,1]), &
matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA)
MatB = MatB + &
matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA)
FAvg = FAvg + F
BMatAvg = BMatAvg + BMat
else
K_eA = K_eA + matmul(transpose(BMat),MatA)
endif
enddo
if (BBarStabilisation) then
FInv = math_inv33(FAvg)
K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + &
(matmul(matmul(transpose(BMatAvg), &
reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + &
K_eB)/real(dimPlex)
else
K_e = K_eA
endif
K_e = K_e + eps*math_identity2nd(cellDof)
K_eVec = reshape(K_e, [cellDof*cellDof])*abs(detJ)
pK_e => K_eVec
call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr)
CHKERRQ(ierr)
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr)
enddo
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! apply boundary conditions
!call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr) MD: linker error
call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_formJacobian
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
use FEM_utilities, only: &
cutBack
use homogenization, only: &
materialpoint_F0, &
materialpoint_F
use FEM_utilities, only: &
utilities_projectBCValues
implicit none
type(tFieldBC), intent(in) :: &
fieldBC
real(pReal), intent(in) :: &
timeinc_old, &
timeinc
logical, intent(in) :: &
guess
PetscInt :: field, face
DM :: dm_local
Vec :: x_local
PetscSection :: section
PetscInt :: bcSize
IS :: bcPoints
PetscErrorCode :: ierr
!--------------------------------------------------------------------------------------------------
! forward last inc
if (guess .and. .not. cutBack) then
ForwardData = .True.
materialpoint_F0 = materialpoint_F
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local
call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0,ierr); CHKERRQ(ierr)
call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector
CHKERRQ(ierr)
call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr)
CHKERRQ(ierr)
call VecAXPY(solution_local,1.0,x_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, &
0.0,fieldBC%componentBC(field)%Value(face),timeinc_old)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
enddo; enddo
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
call VecCopy(solution,solution_rate,ierr); CHKERRQ(ierr)
call VecScale(solution_rate,1.0/timeinc_old,ierr); CHKERRQ(ierr)
endif
call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr)
call VecScale(solution,timeinc,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_forward
!--------------------------------------------------------------------------------------------------
!> @brief reporting
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
err_struct_tolAbs, &
err_struct_tolRel
use IO, only: &
IO_intOut
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: xnorm,snorm,fnorm,divTol
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
!--------------------------------------------------------------------------------------------------
! report
divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*err_struct_tolRel,err_struct_tolAbs)
call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr)
CHKERRQ(ierr)
if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
if (worldrank == 0) then
write(6,'(1/,1x,a,a,i0,a,i0,f0.3)') trim(incInfo), &
' @ Iteration ',PETScIter,' mechanical residual norm = ', &
int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
transpose(P_av)*1.e-6_pReal
flush(6)
endif
end subroutine FEM_mech_converged
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_destroy()
implicit none
PetscErrorCode :: ierr
call VecDestroy(solution,ierr); CHKERRQ(ierr)
call VecDestroy(solution_rate,ierr); CHKERRQ(ierr)
call SNESDestroy(mech_snes,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_destroy
end module FEM_mech

751
src/FEM_utilities.f90 Normal file
View File

@ -0,0 +1,751 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Utilities used by the FEM solver
!--------------------------------------------------------------------------------------------------
module FEM_utilities
#include <petsc/finclude/petsc.h>
use prec, only: pReal, pInt
use PETScdmda
use PETScis
implicit none
private
!--------------------------------------------------------------------------------------------------
!
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer(pInt), public, parameter :: maxFields = 6_pInt
integer(pInt), public :: nActiveFields = 0_pInt
!--------------------------------------------------------------------------------------------------
! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems
real(pReal), public :: wgtDof !< weighting factor 1/Nelems
real(pReal), public :: C_volAvg(3,3,3,3)
!--------------------------------------------------------------------------------------------------
! output data
PetscViewer, public :: resUnit
Vec, public :: coordinatesVec
Vec, allocatable, public :: homogenizationResultsVec(:), &
crystalliteResultsVec(:,:), &
phaseResultsVec(:,:)
!--------------------------------------------------------------------------------------------------
! field labels information
character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical', &
FIELD_THERMAL_label = 'thermal', &
FIELD_DAMAGE_label = 'damage', &
FIELD_SOLUTE_label = 'solute', &
FIELD_MGTWIN_label = 'mgtwin'
enum, bind(c)
enumerator :: FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID
end enum
enum, bind(c)
enumerator :: COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, &
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID
end enum
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical, private :: &
debugGeneral, & !< general debugging of FEM solver
debugRotation, & !< also printing out results in lab frame
debugPETSc !< use some in debug defined options for more verbose PETSc solution
!--------------------------------------------------------------------------------------------------
! derived types
type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true.
logical :: stagConverged = .true.
logical :: regrid = .false.
integer(pInt) :: iterationsNeeded = 0_pInt
end type tSolutionState
type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable :: Value(:)
logical, allocatable :: Mask(:)
end type tComponentBC
type, public :: tFieldBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer(pInt) :: nComponents = 0_pInt
type(tComponentBC), allocatable :: componentBC(:)
end type tFieldBC
type, public :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer(pInt) :: incs = 0_pInt, & !< number of increments
outputfrequency = 1_pInt, & !< frequency of result writes
restartfrequency = 0_pInt, & !< frequency of restart writes
logscale = 0_pInt !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer(pInt), allocatable :: faceID(:)
type(tFieldBC), allocatable :: fieldBC(:)
end type tLoadCase
type, public :: tFEMInterpolation
integer(pInt) :: n
real(pReal), dimension(:,:) , allocatable :: shapeFunc, shapeDerivReal, geomShapeDerivIso
real(pReal), dimension(:,:,:), allocatable :: shapeDerivIso
end type tFEMInterpolation
type, public :: tQuadrature
integer(pInt) :: n
real(pReal), dimension(:) , allocatable :: Weights
real(pReal), dimension(:,:), allocatable :: Points
end type tQuadrature
public :: &
utilities_init, &
utilities_constitutiveResponse, &
utilities_indexBoundaryDofs, &
utilities_projectBCValues, &
utilities_indexActiveSet, &
utilities_destroy, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, &
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID
external :: &
PetscOptionsInsertString, &
PetscObjectSetName, &
DMPlexGetHeightStratum, &
DMGetLabelIdIS, &
DMPlexGetChart, &
DMPlexLabelComplete, &
PetscViewerHDF5Open, &
PetscViewerHDF5PushGroup, &
PetscViewerHDF5PopGroup, &
PetscViewerDestroy
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, sets debug flags
!--------------------------------------------------------------------------------------------------
subroutine utilities_init()
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use DAMASK_interface, only: &
getSolverJobName
use IO, only: &
IO_error, &
IO_warning, &
IO_timeStamp, &
IO_open_file
use numerics, only: &
integrationOrder, &
worldsize, &
worldrank, &
petsc_defaultOptions, &
petsc_options
use debug, only: &
debug_level, &
debug_SPECTRAL, &
debug_LEVELBASIC, &
debug_SPECTRALPETSC, &
debug_SPECTRALROTATION
use debug, only: &
PETSCDEBUG
use math ! must use the whole module for use of FFTW
use mesh, only: &
mesh_NcpElemsGlobal, &
mesh_maxNips, &
geomMesh
use material, only: &
material_homog
implicit none
character(len=1024) :: petsc_optionsPhysics, grainStr
integer(pInt) :: dimPlex
integer(pInt) :: headerID = 205_pInt
PetscInt, dimension(:), pointer :: points
PetscInt, allocatable :: nEntities(:), nOutputCells(:), nOutputNodes(:), mappingCells(:)
PetscInt :: cellStart, cellEnd, cell, ip, dim, ctr, qPt
PetscInt, allocatable :: connectivity(:,:)
Vec :: connectivityVec
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! set debugging parameters
debugGeneral = iand(debug_level(debug_SPECTRAL),debug_LEVELBASIC) /= 0
debugRotation = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION) /= 0
debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0
if(debugPETSc) write(6,'(3(/,a),/)') &
' Initializing PETSc with debug options: ', &
trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.config '
flush(6)
call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
!write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order ' , structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
call PetscViewerHDF5Open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.h5', &
FILE_MODE_WRITE, resUnit, ierr); CHKERRQ(ierr)
call PetscViewerHDF5PushGroup(resUnit, '/', ierr); CHKERRQ(ierr)
call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRQ(ierr)
allocate(nEntities(dimPlex+1), source=0)
allocate(nOutputNodes(worldsize), source = 0)
allocate(nOutputCells(worldsize), source = 0)
do dim = 0, dimPlex
call DMGetStratumSize(geomMesh,'depth',dim,nEntities(dim+1),ierr)
CHKERRQ(ierr)
enddo
select case (integrationOrder)
case(1_pInt)
nOutputNodes(worldrank+1) = nEntities(1)
case(2_pInt)
nOutputNodes(worldrank+1) = sum(nEntities)
case default
nOutputNodes(worldrank+1) = mesh_maxNips*nEntities(dimPlex+1)
end select
nOutputCells(worldrank+1) = count(material_homog > 0_pInt)
call MPI_Allreduce(MPI_IN_PLACE,nOutputNodes,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,nOutputCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (worldrank == 0_pInt) then
open(unit=headerID, file=trim(getSolverJobName())//'.header', &
form='FORMATTED', status='REPLACE')
write(headerID, '(a,i0)') 'dimension : ', dimPlex
write(headerID, '(a,i0)') 'number of nodes : ', sum(nOutputNodes)
write(headerID, '(a,i0)') 'number of cells : ', sum(nOutputCells)
endif
allocate(connectivity(2**dimPlex,nOutputCells(worldrank+1)))
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr)
ctr = 0
select case (integrationOrder)
case(1_pInt)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexGetTransitiveClosure(geomMesh,cell,PETSC_TRUE,points,ierr)
CHKERRQ(ierr)
if (dimPlex == 2) then
connectivity(:,ctr+1) = [points( 9), points(11), points(13), points(13)] - nEntities(dimPlex+1)
ctr = ctr + 1
else
connectivity(:,ctr+1) = [points(23), points(25), points(27), points(27), &
points(29), points(29), points(29), points(29)] - nEntities(dimPlex+1)
ctr = ctr + 1
endif
enddo
case(2_pInt)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexGetTransitiveClosure(geomMesh,cell,PETSC_TRUE,points,ierr)
CHKERRQ(ierr)
if (dimPlex == 2) then
connectivity(:,ctr+1) = [points(9 ), points(3), points(1), points(7)]
connectivity(:,ctr+2) = [points(11), points(5), points(1), points(3)]
connectivity(:,ctr+3) = [points(13), points(7), points(1), points(5)]
ctr = ctr + 3
else
connectivity(:,ctr+1) = [points(23), points(11), points(3), points(15), points(17), points(5), points(1), points(7)]
connectivity(:,ctr+2) = [points(25), points(13), points(3), points(11), points(19), points(9), points(1), points(5)]
connectivity(:,ctr+3) = [points(27), points(15), points(3), points(13), points(21), points(7), points(1), points(9)]
connectivity(:,ctr+4) = [points(29), points(17), points(7), points(21), points(19), points(5), points(1), points(9)]
ctr = ctr + 4_pInt
endif
enddo
case default
do cell = cellStart, cellEnd-1; do ip = 0, mesh_maxNips-1
connectivity(:,ctr+1) = cell*mesh_maxNips + ip
ctr = ctr + 1
enddo; enddo
end select
connectivity = connectivity + sum(nOutputNodes(1:worldrank))
call VecCreateMPI(PETSC_COMM_WORLD,dimPlex*nOutputNodes(worldrank+1),dimPlex*sum(nOutputNodes), &
coordinatesVec,ierr);CHKERRQ(ierr)
call PetscObjectSetName(coordinatesVec, 'NodalCoordinates',ierr)
call VecSetFromOptions(coordinatesVec, ierr); CHKERRQ(ierr)
!allocate(mappingCells(worldsize), source = 0)
!do homog = 1, material_Nhomogenization
! mappingCells = 0_pInt; mappingCells(worldrank+1) = homogOutput(homog)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! call PetscObjectSetName(connectivityVec,'mapping_'//trim(homogenization_name(homog)),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec,results,ierr); CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (material_homog(qPt,cell+1) == homog) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo
!do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
! mappingCells = 0_pInt
! mappingCells(worldrank+1) = crystalliteOutput(cryst,grain)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! call PetscObjectSetName(connectivityVec,'mapping_'// &
! trim(crystallite_name(cryst))//'_'// &
! trim(grainStr),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (homogenization_Ngrains (mesh_element(3,cell+1)) >= grain .and. &
! microstructure_crystallite(mesh_element(4,cell+1)) == cryst) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo; enddo
!do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
! mappingCells = 0_pInt
! mappingCells(worldrank+1) = phaseOutput(phase,grain)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! call PetscObjectSetName(connectivityVec,&
! 'mapping_'//trim(phase_name(phase))//'_'// &
! trim(grainStr),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec, results, ierr)
! CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (material_phase(grain,qPt,cell+1) == phase) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr)
! CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr);CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr);CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr);CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo; enddo
!if (worldrank == 0_pInt) then
! do homog = 1, material_Nhomogenization
! call VecGetSize(homogenizationResultsVec(homog),mappingCells(1),ierr)
! CHKERRQ(ierr)
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of homog_'// &
! trim(homogenization_name(homog))//'_'// &
! 'cells : ', mappingCells(1)
! enddo
! do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
! call VecGetSize(crystalliteResultsVec(cryst,grain),mappingCells(1),ierr)
! CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of cryst_'// &
! trim(crystallite_name(cryst))//'_'// &
! trim(grainStr)//'_'// &
! 'cells : ', mappingCells(1)
! enddo; enddo
! do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
! call VecGetSize(phaseResultsVec(phase,grain),mappingCells(1),ierr)
! CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of phase_'// &
! trim(phase_name(phase))//'_'//trim(grainStr)//'_'// &
! 'cells : ', mappingCells(1)
! enddo; enddo
! close(headerID)
!endif
end subroutine utilities_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates constitutive response
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
use debug, only: &
debug_reset, &
debug_info
use numerics, only: &
worldrank
use math, only: &
math_transpose33, &
math_rotate_forward33, &
math_det33
use FEsolving, only: &
restartWrite
use homogenization, only: &
materialpoint_F0, &
materialpoint_F, &
materialpoint_P, &
materialpoint_dPdF, &
materialpoint_stressAndItsTangent
use mesh, only: &
mesh_NcpElems
implicit none
real(pReal), intent(in) :: timeinc !< loading time
logical, intent(in) :: forwardData !< age results
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
logical :: &
age
integer(pInt) :: &
j
real(pReal) :: defgradDetMin, defgradDetMax, defgradDet
PetscErrorCode :: ierr
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
age = .False.
if (forwardData) then ! aging results
age = .True.
endif
if (cutBack) then ! restore saved variables
age = .False.
endif
call debug_reset()
!--------------------------------------------------------------------------------------------------
! calculate bounds of det(F) and report
if(debugGeneral) then
defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal)
do j = 1_pInt, mesh_NcpElems
defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
end do
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
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
call debug_info()
restartWrite = .false. ! reset restartWrite status
cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,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)
end subroutine utilities_constitutiveResponse
!--------------------------------------------------------------------------------------------------
!> @brief Create index sets of boundary dofs (in local and global numbering)
!--------------------------------------------------------------------------------------------------
subroutine utilities_indexBoundaryDofs(dm_local,nFaceSets,numFields,local2global,section,localIS,globalIS)
implicit none
DM :: dm_local
ISLocalToGlobalMapping :: local2global
PetscSection :: section
PetscInt :: nFaceSets, numFields, nDof
IS, dimension(nFaceSets,numFields) :: localIS, globalIS
PetscInt :: field, faceSet, point, dof, offset
PetscInt :: localSize, storageSize, ISSize
PetscInt, dimension(:) , allocatable :: localIndices
IS :: faceSetIS, BC_IS, dummyIS
PetscInt, dimension(:) , pointer :: pFaceSets, pBCvertex, pBCvertexlc
DMLabel :: BCLabel
PetscErrorCode :: ierr
call DMGetLabel(dm_local,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(dm_local,BCLabel,ierr); CHKERRQ(ierr)
call PetscSectionGetStorageSize(section,storageSize,ierr); CHKERRQ(ierr)
call DMGetLabelIdIS(dm_local,'Face Sets',faceSetIS,ierr); CHKERRQ(ierr)
call ISGetIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
allocate(localIndices (storageSize))
do faceSet = 1, nFaceSets
call DMGetStratumSize(dm_local,'Face Sets',pFaceSets(faceSet),ISSize,ierr)
CHKERRQ(ierr)
call DMGetStratumIS(dm_local,'Face Sets',pFaceSets(faceSet),BC_IS,ierr)
CHKERRQ(ierr)
if (ISSize > 0) call ISGetIndicesF90(BC_IS,pBCvertex,ierr)
do field = 1, numFields
localSize = 0
do point = 1, ISSize
call PetscSectionGetFieldDof(section,pBCvertex(point),field-1,nDof,ierr)
CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,pBCvertex(point),field-1,offset,ierr)
CHKERRQ(ierr)
do dof = 1, nDof
localSize = localSize + 1
localIndices(localSize) = offset + dof - 1
enddo
enddo
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
localIS(faceSet,field),ierr)
CHKERRQ(ierr)
call ISLocalToGlobalMappingApplyIS(local2global,localIS(faceSet,field), &
globalIS(faceSet,field),ierr)
CHKERRQ(ierr)
enddo
if (ISSize > 0) call ISRestoreIndicesF90(BC_IS,pBCvertex,ierr)
call ISDestroy(BC_IS,ierr); CHKERRQ(ierr)
enddo
call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
call ISDestroy(faceSetIS,ierr); CHKERRQ(ierr)
do faceSet = 1, nFaceSets; do field = 1, numFields
call ISGetSize(globalIS(faceSet,field),ISSize,ierr); CHKERRQ(ierr)
if (ISSize > 0) then
call ISGetIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
call ISGetIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
endif
localSize = 0
do point = 1, ISSize
if (pBCvertex(point) >= 0) then
localSize = localSize + 1
localIndices(localSize) = pBCvertexlc(point)
endif
enddo
if (ISSize > 0) then
call ISRestoreIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
call ISRestoreIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
endif
call ISDestroy(globalIS(faceSet,field),ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
globalIS(faceSet,field),ierr)
CHKERRQ(ierr)
if (ISSize > 0) then
call ISDuplicate(localIS(faceSet,field),dummyIS,ierr); CHKERRQ(ierr)
call ISDestroy(localIS(faceSet,field),ierr); CHKERRQ(ierr)
call ISDifference(dummyIS,globalIS(faceSet,field),localIS(faceSet,field),ierr)
CHKERRQ(ierr)
call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
endif
enddo; enddo
deallocate(localIndices)
end subroutine utilities_indexBoundaryDofs
!--------------------------------------------------------------------------------------------------
!> @brief Project BC values to local vector
!--------------------------------------------------------------------------------------------------
subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc)
implicit none
Vec :: localVec
PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset
PetscSection :: section
IS :: bcPointsIS
PetscInt, pointer :: bcPoints(:)
PetscScalar, pointer :: localArray(:)
PetscScalar :: BCValue,BCDotValue,timeinc
PetscErrorCode :: ierr
call PetscSectionGetFieldComponents(section,field,numComp,ierr); CHKERRQ(ierr)
call ISGetSize(bcPointsIS,nBcPoints,ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISGetIndicesF90(bcPointsIS,bcPoints,ierr)
call VecGetArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
do point = 1, nBcPoints
call PetscSectionGetFieldDof(section,bcPoints(point),field,numDof,ierr)
CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,ierr)
CHKERRQ(ierr)
do dof = offset+comp+1, offset+numDof, numComp
localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
enddo
enddo
call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr)
call VecAssemblyEnd (localVec, ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr)
end subroutine utilities_projectBCValues
!--------------------------------------------------------------------------------------------------
!> @brief Create index sets of boundary dofs (in local and global numbering)
!--------------------------------------------------------------------------------------------------
subroutine utilities_indexActiveSet(field,section,x_local,f_local,localIS,globalIS)
use mesh, only: &
geomMesh
implicit none
ISLocalToGlobalMapping :: local2global
PetscSection :: section
Vec :: x_local, f_local
PetscInt :: field
IS :: localIS, globalIS, dummyIS
PetscScalar, dimension(:) , pointer :: x_scal, f_scal
PetscInt :: ISSize
PetscInt :: chart, chartStart, chartEnd, nDof, dof, offset
PetscInt :: localSize
PetscInt, dimension(:) , allocatable :: localIndices
PetscInt, dimension(:) , pointer :: pBCvertex, pBCvertexlc
PetscErrorCode :: ierr
call DMGetLocalToGlobalMapping(geomMesh,local2global,ierr)
CHKERRQ(ierr)
call DMPlexGetChart(geomMesh,chartStart,chartEnd,ierr)
CHKERRQ(ierr)
call VecGetArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
call VecGetArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
localSize = 0
do chart = chartStart, chartEnd-1
call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
do dof = offset+1, offset+nDof
if (((x_scal(dof) < 1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) localSize = localSize + 1
enddo
enddo
allocate(localIndices(localSize))
localSize = 0
do chart = chartStart, chartEnd-1
call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
do dof = offset+1, offset+nDof
if (((x_scal(dof) < 1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) then
localSize = localSize + 1
localIndices(localSize) = dof-1
endif
enddo
enddo
call VecRestoreArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
call VecRestoreArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,localIS,ierr)
CHKERRQ(ierr)
call ISLocalToGlobalMappingApplyIS(local2global,localIS,globalIS,ierr)
CHKERRQ(ierr)
call ISGetSize(globalIS,ISSize,ierr); CHKERRQ(ierr)
if (ISSize > 0) then
call ISGetIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
call ISGetIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
endif
localSize = 0
do chart = 1, ISSize
if (pBCvertex(chart) >= 0) then
localSize = localSize + 1
localIndices(localSize) = pBCvertexlc(chart)
endif
enddo
if (ISSize > 0) then
call ISRestoreIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
call ISRestoreIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
endif
call ISDestroy(globalIS,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,globalIS,ierr)
CHKERRQ(ierr)
if (ISSize > 0) then
call ISDuplicate(localIS,dummyIS,ierr); CHKERRQ(ierr)
call ISDestroy(localIS,ierr); CHKERRQ(ierr)
call ISDifference(dummyIS,globalIS,localIS,ierr)
CHKERRQ(ierr)
call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
endif
deallocate(localIndices)
end subroutine utilities_indexActiveSet
!--------------------------------------------------------------------------------------------------
!> @brief cleans up
!--------------------------------------------------------------------------------------------------
subroutine utilities_destroy()
!use material, only: &
! homogenization_Ngrains
!implicit none
!PetscInt :: homog, cryst, grain, phase
!PetscErrorCode :: ierr
!call PetscViewerHDF5PopGroup(resUnit, ierr); CHKERRQ(ierr)
!call VecDestroy(coordinatesVec,ierr); CHKERRQ(ierr)
!do homog = 1, material_Nhomogenization
! call VecDestroy(homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
! do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_Ngrains(homog)
! call VecDestroy(crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
! enddo; enddo
! do phase = 1, material_Nphase; do grain = 1, homogenization_Ngrains(homog)
! call VecDestroy(phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
! enddo; enddo
!enddo
!call PetscViewerDestroy(resUnit, ierr); CHKERRQ(ierr)
end subroutine utilities_destroy
end module FEM_utilities

350
src/FEM_zoo.f90 Normal file
View File

@ -0,0 +1,350 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Interpolation data used by the FEM solver
!--------------------------------------------------------------------------------------------------
module FEM_Zoo
use prec, only: pReal, pInt, group_float
implicit none
private
integer(pInt), parameter, public:: &
maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary)
real(pReal), dimension(2,3), private, protected :: &
triangle = reshape([-1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal], shape=[2,3])
real(pReal), dimension(3,4), private, protected :: &
tetrahedron = reshape([-1.0_pReal, -1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal, -1.0_pReal, &
-1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4])
integer(pInt), dimension(3,maxOrder), public, protected :: &
FEM_Zoo_nQuadrature !< number of quadrature points for a given spatial dimension(1-3) and interpolation order(1-maxOrder)
type(group_float), dimension(3,maxOrder), public, protected :: &
FEM_Zoo_QuadratureWeights, & !< quadrature weights for each quadrature rule
FEM_Zoo_QuadraturePoints !< quadrature point coordinates (in simplical system) for each quadrature rule
public :: &
FEM_Zoo_init
contains
!--------------------------------------------------------------------------------------------------
!> @brief initializes FEM interpolation data
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_timeStamp
implicit none
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'
write(6,'(a)') ' $Id: FEM_Zoo.f90 4354 2015-08-04 15:04:53Z MPIE\p.shanthraj $'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! 2D linear
FEM_Zoo_nQuadrature(2,1) = 1
allocate(FEM_Zoo_QuadratureWeights(2,1)%p(1))
allocate(FEM_Zoo_QuadraturePoints (2,1)%p(2))
FEM_Zoo_QuadratureWeights(2,1)%p(1) = 1.0_pReal
call FEM_Zoo_permutationStar3([1.0_pReal/3.0_pReal], &
FEM_Zoo_QuadraturePoints(2,1)%p(1:2))
!--------------------------------------------------------------------------------------------------
! 2D quadratic
FEM_Zoo_nQuadrature(2,2) = 3
allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3))
allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6))
FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
call FEM_Zoo_permutationStar21([1.0_pReal/6.0_pReal], &
FEM_Zoo_QuadraturePoints(2,2)%p(1:6))
!--------------------------------------------------------------------------------------------------
! 2D cubic
FEM_Zoo_nQuadrature(2,3) = 6
allocate(FEM_Zoo_QuadratureWeights(2,3)%p(6 ))
allocate(FEM_Zoo_QuadraturePoints (2,3)%p(12))
FEM_Zoo_QuadratureWeights(2,3)%p(1:3) = 0.22338158967801146570_pReal
call FEM_Zoo_permutationStar21([0.44594849091596488632_pReal], &
FEM_Zoo_QuadraturePoints(2,3)%p(1:6))
FEM_Zoo_QuadratureWeights(2,3)%p(4:6) = 0.10995174365532186764_pReal
call FEM_Zoo_permutationStar21([0.091576213509770743460_pReal], &
FEM_Zoo_QuadraturePoints(2,3)%p(7:12))
!--------------------------------------------------------------------------------------------------
! 2D quartic
FEM_Zoo_nQuadrature(2,4) = 12
allocate(FEM_Zoo_QuadratureWeights(2,4)%p(12))
allocate(FEM_Zoo_QuadraturePoints (2,4)%p(24))
FEM_Zoo_QuadratureWeights(2,4)%p(1:3) = 0.11678627572638_pReal
call FEM_Zoo_permutationStar21([0.24928674517091_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(1:6))
FEM_Zoo_QuadratureWeights(2,4)%p(4:6) = 0.05084490637021_pReal
call FEM_Zoo_permutationStar21([0.06308901449150_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(7:12))
FEM_Zoo_QuadratureWeights(2,4)%p(7:12) = 0.08285107561837_pReal
call FEM_Zoo_permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(13:24))
!--------------------------------------------------------------------------------------------------
! 2D order 5
FEM_Zoo_nQuadrature(2,5) = 16
allocate(FEM_Zoo_QuadratureWeights(2,5)%p(16))
allocate(FEM_Zoo_QuadraturePoints (2,5)%p(32))
FEM_Zoo_QuadratureWeights(2,5)%p(1 ) = 0.14431560767779_pReal
call FEM_Zoo_permutationStar3([0.33333333333333_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(1:2))
FEM_Zoo_QuadratureWeights(2,5)%p(2:4) = 0.09509163426728_pReal
call FEM_Zoo_permutationStar21([0.45929258829272_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(3:8))
FEM_Zoo_QuadratureWeights(2,5)%p(5:7) = 0.10321737053472_pReal
call FEM_Zoo_permutationStar21([0.17056930775176_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(9:14))
FEM_Zoo_QuadratureWeights(2,5)%p(8:10) = 0.03245849762320_pReal
call FEM_Zoo_permutationStar21([0.05054722831703_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(15:20))
FEM_Zoo_QuadratureWeights(2,5)%p(11:16) = 0.02723031417443_pReal
call FEM_Zoo_permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(21:32))
!--------------------------------------------------------------------------------------------------
! 3D linear
FEM_Zoo_nQuadrature(3,1) = 1
allocate(FEM_Zoo_QuadratureWeights(3,1)%p(1))
allocate(FEM_Zoo_QuadraturePoints (3,1)%p(3))
FEM_Zoo_QuadratureWeights(3,1)%p(1) = 1.0_pReal
call FEM_Zoo_permutationStar4([0.25_pReal], &
FEM_Zoo_QuadraturePoints(3,1)%p(1:3))
!--------------------------------------------------------------------------------------------------
! 3D quadratic
FEM_Zoo_nQuadrature(3,2) = 4
allocate(FEM_Zoo_QuadratureWeights(3,2)%p(4 ))
allocate(FEM_Zoo_QuadraturePoints (3,2)%p(12))
FEM_Zoo_QuadratureWeights(3,2)%p(1:4) = 0.25_pReal
call FEM_Zoo_permutationStar31([0.13819660112501051518_pReal], &
FEM_Zoo_QuadraturePoints(3,2)%p(1:12))
!--------------------------------------------------------------------------------------------------
! 3D cubic
FEM_Zoo_nQuadrature(3,3) = 14
allocate(FEM_Zoo_QuadratureWeights(3,3)%p(14))
allocate(FEM_Zoo_QuadraturePoints (3,3)%p(42))
FEM_Zoo_QuadratureWeights(3,3)%p(1:4) = 0.073493043116361949544_pReal
call FEM_Zoo_permutationStar31([0.092735250310891226402_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(1:12))
FEM_Zoo_QuadratureWeights(3,3)%p(5:8) = 0.11268792571801585080_pReal
call FEM_Zoo_permutationStar31([0.31088591926330060980_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(13:24))
FEM_Zoo_QuadratureWeights(3,3)%p(9:14) = 0.042546020777081466438_pReal
call FEM_Zoo_permutationStar22([0.045503704125649649492_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(25:42))
!--------------------------------------------------------------------------------------------------
! 3D quartic
FEM_Zoo_nQuadrature(3,4) = 35
allocate(FEM_Zoo_QuadratureWeights(3,4)%p(35))
allocate(FEM_Zoo_QuadraturePoints (3,4)%p(105))
FEM_Zoo_QuadratureWeights(3,4)%p(1:4) = 0.0021900463965388_pReal
call FEM_Zoo_permutationStar31([0.0267367755543735_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(1:12))
FEM_Zoo_QuadratureWeights(3,4)%p(5:16) = 0.0143395670177665_pReal
call FEM_Zoo_permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(13:48))
FEM_Zoo_QuadratureWeights(3,4)%p(17:22) = 0.0250305395686746_pReal
call FEM_Zoo_permutationStar22([0.4547545999844830_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(49:66))
FEM_Zoo_QuadratureWeights(3,4)%p(23:34) = 0.0479839333057554_pReal
call FEM_Zoo_permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(67:102))
FEM_Zoo_QuadratureWeights(3,4)%p(35) = 0.0931745731195340_pReal
call FEM_Zoo_permutationStar4([0.25_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(103:105))
!--------------------------------------------------------------------------------------------------
! 3D quintic
FEM_Zoo_nQuadrature(3,5) = 56
allocate(FEM_Zoo_QuadratureWeights(3,5)%p(56))
allocate(FEM_Zoo_QuadraturePoints (3,5)%p(168))
FEM_Zoo_QuadratureWeights(3,5)%p(1:4) = 0.0010373112336140_pReal
call FEM_Zoo_permutationStar31([0.0149520651530592_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(1:12))
FEM_Zoo_QuadratureWeights(3,5)%p(5:16) = 0.0096016645399480_pReal
call FEM_Zoo_permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(13:48))
FEM_Zoo_QuadratureWeights(3,5)%p(17:28) = 0.0164493976798232_pReal
call FEM_Zoo_permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(49:84))
FEM_Zoo_QuadratureWeights(3,5)%p(29:40) = 0.0153747766513310_pReal
call FEM_Zoo_permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(85:120))
FEM_Zoo_QuadratureWeights(3,5)%p(41:52) = 0.0293520118375230_pReal
call FEM_Zoo_permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(121:156))
FEM_Zoo_QuadratureWeights(3,5)%p(53:56) = 0.0366291366405108_pReal
call FEM_Zoo_permutationStar31([0.1344783347929940_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(157:168))
end subroutine FEM_Zoo_init
!--------------------------------------------------------------------------------------------------
!> @brief star 3 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar3(point,qPt)
implicit none
real(pReal) :: point(1), qPt(2,1), temp(3,1)
temp(:,1) = [point(1), point(1), point(1)]
qPt = matmul(triangle, temp)
end subroutine FEM_Zoo_permutationStar3
!--------------------------------------------------------------------------------------------------
!> @brief star 21 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar21(point,qPt)
implicit none
real(pReal) :: point(1), qPt(2,3), temp(3,3)
temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)]
temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)]
temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)]
qPt = matmul(triangle, temp)
end subroutine FEM_Zoo_permutationStar21
!--------------------------------------------------------------------------------------------------
!> @brief star 111 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar111(point,qPt)
implicit none
real(pReal) :: point(2), qPt(2,6), temp(3,6)
temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)]
temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)]
temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)]
temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)]
temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)]
qPt = matmul(triangle, temp)
end subroutine FEM_Zoo_permutationStar111
!--------------------------------------------------------------------------------------------------
!> @brief star 4 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar4(point,qPt)
implicit none
real(pReal) :: point(1), qPt(3,1), temp(4,1)
temp(:,1) = [point(1), point(1), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar4
!--------------------------------------------------------------------------------------------------
!> @brief star 31 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar31(point,qPt)
implicit none
real(pReal) :: point(1), qPt(3,4), temp(4,4)
temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)]
temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)]
temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)]
temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar31
!--------------------------------------------------------------------------------------------------
!> @brief star 22 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar22(point,qPt)
implicit none
real(pReal) :: point(1), qPt(3,6), temp(4,6)
temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)]
temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)]
temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)]
temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)]
temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)]
temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar22
!--------------------------------------------------------------------------------------------------
!> @brief star 211 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar211(point,qPt)
implicit none
real(pReal) :: point(2), qPt(3,12), temp(4,12)
temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)]
temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)]
temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)]
temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)]
temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)]
temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)]
temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar211
!--------------------------------------------------------------------------------------------------
!> @brief star 1111 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar1111(point,qPt)
implicit none
real(pReal) :: point(3), qPt(3,24), temp(4,24)
temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)]
temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)]
temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)]
temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)]
temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)]
temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)]
temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)]
temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)]
temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)]
temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)]
temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)]
temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar1111
end module FEM_Zoo

9
src/FEsolving.f90
View File

@ -81,20 +81,13 @@ subroutine FE_init
modelName = getSolverJobName() modelName = getSolverJobName()
#if defined(Spectral) || defined(FEM) #if defined(Spectral) || defined(FEM)
restartInc = interface_RestartInc
#ifdef Spectral
restartInc = spectralRestartInc
#endif
#ifdef FEM
restartInc = FEMRestartInc
#endif
if(restartInc < 0_pInt) then if(restartInc < 0_pInt) then
call IO_warning(warning_ID=34_pInt) call IO_warning(warning_ID=34_pInt)
restartInc = 0_pInt restartInc = 0_pInt
endif endif
restartRead = restartInc > 0_pInt ! only read in if "true" restart requested restartRead = restartInc > 0_pInt ! only read in if "true" restart requested
#else #else
call IO_open_inputFile(FILEUNIT,modelName) call IO_open_inputFile(FILEUNIT,modelName)
rewind(FILEUNIT) rewind(FILEUNIT)

301
src/IO.f90
View File

@ -22,6 +22,7 @@ module IO
public :: & public :: &
IO_init, & IO_init, &
IO_read, & IO_read, &
IO_recursiveRead, &
IO_checkAndRewind, & IO_checkAndRewind, &
IO_open_file_stat, & IO_open_file_stat, &
IO_open_jobFile_stat, & IO_open_jobFile_stat, &
@ -35,10 +36,6 @@ module IO
IO_hybridIA, & IO_hybridIA, &
IO_isBlank, & IO_isBlank, &
IO_getTag, & IO_getTag, &
IO_countSections, &
IO_countTagInPart, &
IO_spotTagInPart, &
IO_globalTagInPart, &
IO_stringPos, & IO_stringPos, &
IO_stringValue, & IO_stringValue, &
IO_fixedStringValue ,& IO_fixedStringValue ,&
@ -100,6 +97,7 @@ end subroutine IO_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief recursively reads a line from a text file. !> @brief recursively reads a line from a text file.
!! Recursion is triggered by "{path/to/inputfile}" in a line !! Recursion is triggered by "{path/to/inputfile}" in a line
!> @details unstable and buggy
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
recursive function IO_read(fileUnit,reset) result(line) recursive function IO_read(fileUnit,reset) result(line)
@ -151,7 +149,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),action='read') ! open included file open(newunit=unitOn(stack),iostat=myStat,file=pathOn(stack),action='read',status='old',position='rewind') ! 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)
@ -170,6 +168,80 @@ recursive function IO_read(fileUnit,reset) result(line)
end function IO_read end function IO_read
!--------------------------------------------------------------------------------------------------
!> @brief recursively reads a text file.
!! Recursion is triggered by "{path/to/inputfile}" in a line
!--------------------------------------------------------------------------------------------------
recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
implicit none
character(len=*), intent(in) :: fileName
integer(pInt), intent(in), optional :: cnt !< recursion counter
character(len=256), dimension(:), allocatable :: fileContent !< file content, separated per lines
character(len=256), dimension(:), allocatable :: includedContent
character(len=256) :: line
character(len=256), parameter :: dummy = 'https://damask.mpie.de' !< to fill up remaining array
character(len=:), allocatable :: rawData
integer(pInt) :: &
fileLength, &
fileUnit, &
startPos, endPos, &
myTotalLines, & !< # lines read from file without include statements
includedLines, & !< # lines included from other file(s)
missingLines, & !< # lines missing from current file
l,i, &
myStat
if (merge(cnt,0_pInt,present(cnt))>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName))
!--------------------------------------------------------------------------------------------------
! read data as stream
inquire(file = fileName, size=fileLength)
open(newunit=fileUnit, file=fileName, access='stream',&
status='old', position='rewind', action='read',iostat=myStat)
if(myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=trim(fileName))
allocate(character(len=fileLength)::rawData)
read(fileUnit) rawData
close(fileUnit)
!--------------------------------------------------------------------------------------------------
! count lines to allocate string array
myTotalLines = 0_pInt
do l=1_pInt, len(rawData)
if (rawData(l:l) == new_line('')) myTotalLines = myTotalLines+1
enddo
allocate(fileContent(myTotalLines))
!--------------------------------------------------------------------------------------------------
! split raw data at end of line and handle includes
startPos = 1_pInt
endPos = 0_pInt
includedLines=0_pInt
l=0_pInt
do while (startPos <= len(rawData))
l = l + 1_pInt
endPos = endPos + scan(rawData(startPos:),new_line(''))
if(endPos - startPos >256) call IO_error(107_pInt,ext_msg=trim(fileName))
line = rawData(startPos:endPos-1_pInt)
startPos = endPos + 1_pInt
recursion: if(scan(trim(line),'{') < scan(trim(line),'}')) then
myTotalLines = myTotalLines - 1_pInt
includedContent = IO_recursiveRead(trim(line(scan(line,'{')+1_pInt:scan(line,'}')-1_pInt)), &
merge(cnt,1_pInt,present(cnt))) ! to track recursion depth
includedLines = includedLines + size(includedContent)
missingLines = myTotalLines + includedLines - size(fileContent(1:l-1)) -size(includedContent)
fileContent = [ fileContent(1:l-1_pInt), includedContent, [(dummy,i=1,missingLines)] ] ! add content and grow array
l = l - 1_pInt + size(includedContent)
else recursion
fileContent(l) = line
endif recursion
enddo
end function IO_recursiveRead
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief checks if unit is opened for reading, if true rewinds. Otherwise stops with !> @brief checks if unit is opened for reading, if true rewinds. Otherwise stops with
@ -178,7 +250,7 @@ end function IO_read
subroutine IO_checkAndRewind(fileUnit) subroutine IO_checkAndRewind(fileUnit)
implicit none implicit none
integer(pInt), intent(in) :: fileUnit !< file unit integer(pInt), intent(in) :: fileUnit !< file unit
logical :: fileOpened logical :: fileOpened
character(len=15) :: fileRead character(len=15) :: fileRead
@ -203,7 +275,7 @@ subroutine IO_open_file(fileUnit,path)
integer(pInt) :: myStat integer(pInt) :: myStat
open(fileUnit,status='old',iostat=myStat,file=path) open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
end subroutine IO_open_file end subroutine IO_open_file
@ -222,7 +294,8 @@ logical function IO_open_file_stat(fileUnit,path)
integer(pInt) :: myStat integer(pInt) :: myStat
open(fileUnit,status='old',iostat=myStat,file=path) open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) close(fileUnit)
IO_open_file_stat = (myStat == 0_pInt) IO_open_file_stat = (myStat == 0_pInt)
end function IO_open_file_stat end function IO_open_file_stat
@ -246,7 +319,7 @@ subroutine IO_open_jobFile(fileUnit,ext)
character(len=1024) :: path character(len=1024) :: path
path = trim(getSolverJobName())//'.'//ext path = trim(getSolverJobName())//'.'//ext
open(fileUnit,status='old',iostat=myStat,file=path) open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
end subroutine IO_open_jobFile end subroutine IO_open_jobFile
@ -270,7 +343,8 @@ logical function IO_open_jobFile_stat(fileUnit,ext)
character(len=1024) :: path character(len=1024) :: path
path = trim(getSolverJobName())//'.'//ext path = trim(getSolverJobName())//'.'//ext
open(fileUnit,status='old',iostat=myStat,file=path) open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) close(fileUnit)
IO_open_jobFile_stat = (myStat == 0_pInt) IO_open_jobFile_stat = (myStat == 0_pInt)
end function IO_open_JobFile_stat end function IO_open_JobFile_stat
@ -296,11 +370,11 @@ subroutine IO_open_inputFile(fileUnit,modelName)
fileType = 1_pInt ! assume .pes fileType = 1_pInt ! assume .pes
path = trim(modelName)//inputFileExtension(fileType) ! attempt .pes, if it exists: it should be used path = trim(modelName)//inputFileExtension(fileType) ! attempt .pes, if it exists: it should be used
open(fileUnit+1,status='old',iostat=myStat,file=path) open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
if(myStat /= 0_pInt) then ! if .pes does not work / exist; use conventional extension, i.e.".inp" if(myStat /= 0_pInt) then ! if .pes does not work / exist; use conventional extension, i.e.".inp"
fileType = 2_pInt fileType = 2_pInt
path = trim(modelName)//inputFileExtension(fileType) path = trim(modelName)//inputFileExtension(fileType)
open(fileUnit+1,status='old',iostat=myStat,file=path) open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
endif endif
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
@ -335,7 +409,7 @@ subroutine IO_open_logFile(fileUnit)
character(len=1024) :: path character(len=1024) :: path
path = trim(getSolverJobName())//LogFileExtension path = trim(getSolverJobName())//LogFileExtension
open(fileUnit,status='old',iostat=myStat,file=path) open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path) if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
end subroutine IO_open_logFile end subroutine IO_open_logFile
@ -755,16 +829,22 @@ pure function IO_getTag(string,openChar,closeChar)
character(len=*), intent(in) :: string !< string to check for tag character(len=*), intent(in) :: string !< string to check for tag
character(len=len_trim(string)) :: IO_getTag character(len=len_trim(string)) :: IO_getTag
character(len=*), intent(in) :: openChar, & !< indicates beginning of tag character, intent(in) :: openChar, & !< indicates beginning of tag
closeChar !< indicates end of tag closeChar !< indicates end of tag
character(len=*), parameter :: SEP=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces character(len=*), parameter :: SEP=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces
integer :: left,right ! no pInt integer :: left,right ! no pInt
IO_getTag = '' IO_getTag = ''
left = scan(string,openChar)
right = scan(string,closeChar)
if (openChar /= closeChar) then
left = scan(string,openChar)
right = scan(string,closeChar)
else
left = scan(string,openChar)
right = left + merge(scan(string(left+1:),openChar),0_pInt,len(string) > left)
endif
if (left == verify(string,SEP) .and. right > left) & ! openChar is first and closeChar occurs if (left == verify(string,SEP) .and. right > left) & ! openChar is first and closeChar occurs
IO_getTag = string(left+1:right-1) IO_getTag = string(left+1:right-1)
@ -772,173 +852,6 @@ pure function IO_getTag(string,openChar,closeChar)
end function IO_getTag end function IO_getTag
!--------------------------------------------------------------------------------------------------
!> @brief count number of [sections] in <part> for given file handle
!--------------------------------------------------------------------------------------------------
integer(pInt) function IO_countSections(fileUnit,part)
implicit none
integer(pInt), intent(in) :: fileUnit !< file handle
character(len=*), intent(in) :: part !< part name in which sections are counted
character(len=65536) :: line
line = ''
IO_countSections = 0_pInt
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part) ! search for part
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier
IO_countSections = IO_countSections + 1_pInt
enddo
end function IO_countSections
!--------------------------------------------------------------------------------------------------
!> @brief returns array of tag counts within <part> for at most N [sections]
!--------------------------------------------------------------------------------------------------
function IO_countTagInPart(fileUnit,part,tag,Nsections)
implicit none
integer(pInt), intent(in) :: Nsections !< maximum number of sections in which tag is searched for
integer(pInt), dimension(Nsections) :: IO_countTagInPart
integer(pInt), intent(in) :: fileUnit !< file handle
character(len=*),intent(in) :: part, & !< part in which tag is searched for
tag !< tag to search for
integer(pInt), dimension(Nsections) :: counter
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: section
character(len=65536) :: line
line = ''
counter = 0_pInt
section = 0_pInt
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part) ! search for part
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') section = section + 1_pInt ! found [section] identifier
if (section > 0) then
chunkPos = IO_stringPos(line)
if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) & ! match
counter(section) = counter(section) + 1_pInt
endif
enddo
IO_countTagInPart = counter
end function IO_countTagInPart
!--------------------------------------------------------------------------------------------------
!> @brief returns array of tag presence within <part> for at most N [sections]
!--------------------------------------------------------------------------------------------------
function IO_spotTagInPart(fileUnit,part,tag,Nsections)
implicit none
integer(pInt), intent(in) :: Nsections !< maximum number of sections in which tag is searched for
logical, dimension(Nsections) :: IO_spotTagInPart
integer(pInt), intent(in) :: fileUnit !< file handle
character(len=*),intent(in) :: part, & !< part in which tag is searched for
tag !< tag to search for
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: section
character(len=65536) :: line
IO_spotTagInPart = .false. ! assume to nowhere spot tag
section = 0_pInt
line = ''
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part) ! search for part
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
foundNextPart: if (IO_getTag(line,'<','>') /= '') then
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif foundNextPart
if (IO_getTag(line,'[',']') /= '') section = section + 1_pInt ! found [section] identifier
if (section > 0_pInt) then
chunkPos = IO_stringPos(line)
if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) & ! match
IO_spotTagInPart(section) = .true.
endif
enddo
end function IO_spotTagInPart
!--------------------------------------------------------------------------------------------------
!> @brief return logical whether tag is present within <part> before any [sections]
!--------------------------------------------------------------------------------------------------
logical function IO_globalTagInPart(fileUnit,part,tag)
implicit none
integer(pInt), intent(in) :: fileUnit !< file handle
character(len=*),intent(in) :: part, & !< part in which tag is searched for
tag !< tag to search for
integer(pInt), allocatable, dimension(:) :: chunkPos
character(len=65536) :: line
IO_globalTagInPart = .false. ! assume to nowhere spot tag
line =''
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part) ! search for part
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
foundNextPart: if (IO_getTag(line,'<','>') /= '') then
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif foundNextPart
foundFirstSection: if (IO_getTag(line,'[',']') /= '') then
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif foundFirstSection
chunkPos = IO_stringPos(line)
match: if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) then
IO_globalTagInPart = .true.
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif match
enddo
end function IO_globalTagInPart
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief locates all space-separated chunks in given string and returns array containing number !> @brief locates all space-separated chunks in given string and returns array containing number
!! them and the left/right position to be used by IO_xxxVal !! them and the left/right position to be used by IO_xxxVal
@ -1450,12 +1363,16 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
pure function IO_intOut(intToPrint) pure function IO_intOut(intToPrint)
implicit none implicit none
character(len=19) :: N_Digits ! maximum digits for 64 bit integer
character(len=40) :: IO_intOut
integer(pInt), intent(in) :: intToPrint integer(pInt), intent(in) :: intToPrint
character(len=41) :: IO_intOut
integer(pInt) :: N_digits
character(len=19) :: width ! maximum digits for 64 bit integer
character(len=20) :: min_width ! longer for negative values
write(N_Digits, '(I19.19)') 1_pInt + int(log10(real(intToPrint)),pInt) N_digits = 1_pInt + int(log10(real(max(abs(intToPrint),1_pInt))),pInt)
IO_intOut = 'I'//trim(N_Digits)//'.'//trim(N_Digits) write(width, '(I19.19)') N_digits
write(min_width, '(I20.20)') N_digits + merge(1_pInt,0_pInt,intToPrint < 0_pInt)
IO_intOut = 'I'//trim(min_width)//'.'//trim(width)
end function IO_intOut end function IO_intOut
@ -1513,6 +1430,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'unknown output:' msg = 'unknown output:'
case (106_pInt) case (106_pInt)
msg = 'working directory does not exist:' msg = 'working directory does not exist:'
case (107_pInt)
msg = 'line length exceeds limit of 256'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! lattice error messages ! lattice error messages
@ -1558,6 +1477,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'illegal texture transformation specified' msg = 'illegal texture transformation specified'
case (160_pInt) case (160_pInt)
msg = 'no entries in config part' msg = 'no entries in config part'
case (161_pInt)
msg = 'config part found twice'
case (165_pInt) case (165_pInt)
msg = 'homogenization configuration' msg = 'homogenization configuration'
case (170_pInt) case (170_pInt)
@ -1655,7 +1576,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
case (845_pInt) case (845_pInt)
msg = 'incomplete information in spectral mesh header' msg = 'incomplete information in spectral mesh header'
case (846_pInt) case (846_pInt)
msg = 'not a rotation defined for loadcase rotation' msg = 'rotation for load case rotation ill-defined (R:RT != I)'
case (847_pInt) case (847_pInt)
msg = 'update of gamma operator not possible when pre-calculated' msg = 'update of gamma operator not possible when pre-calculated'
case (880_pInt) case (880_pInt)

327
src/config.f90
View File

@ -20,12 +20,17 @@ module config
type, public :: tPartitionedStringList type, public :: tPartitionedStringList
type(tPartitionedString) :: string type(tPartitionedString) :: string
type(tPartitionedStringList), pointer :: next => null() type(tPartitionedStringList), pointer :: next => null()
contains contains
procedure :: add => add procedure :: add => add
procedure :: show => show procedure :: show => show
procedure :: free => free procedure :: free => free
! currently, a finalize is needed for all shapes of tPartitionedStringList.
! with Fortran 2015, we can define one recursive elemental function
! https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/543326
final :: finalize, &
finalizeArray
procedure :: keyExists => keyExists procedure :: keyExists => keyExists
procedure :: countKeys => countKeys procedure :: countKeys => countKeys
@ -37,11 +42,10 @@ module config
procedure :: getInts => getInts procedure :: getInts => getInts
procedure :: getStrings => getStrings procedure :: getStrings => getStrings
end type tPartitionedStringList end type tPartitionedStringList
type(tPartitionedStringList), public :: emptyList type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: &
type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: & ! QUESTION: rename to config_XXX?
config_phase, & config_phase, &
config_microstructure, & config_microstructure, &
config_homogenization, & config_homogenization, &
@ -76,7 +80,6 @@ module config
MATERIAL_configFile = 'material.config', & !< generic name for material configuration file MATERIAL_configFile = 'material.config', & !< generic name for material configuration file
MATERIAL_localFileExt = 'materialConfig' !< extension of solver job name depending material configuration file MATERIAL_localFileExt = 'materialConfig' !< extension of solver job name depending material configuration file
public :: & public :: &
config_init, & config_init, &
config_deallocate config_deallocate
@ -92,12 +95,14 @@ subroutine config_init()
compiler_version, & compiler_version, &
compiler_options compiler_options
#endif #endif
use prec, only: &
pStringLen
use DAMASK_interface, only: &
getSolverJobName
use IO, only: & use IO, only: &
IO_error, & IO_error, &
IO_open_file, &
IO_read, &
IO_lc, & IO_lc, &
IO_open_jobFile_stat, & IO_recursiveRead, &
IO_getTag, & IO_getTag, &
IO_timeStamp, & IO_timeStamp, &
IO_EOF IO_EOF
@ -107,12 +112,13 @@ subroutine config_init()
debug_levelBasic debug_levelBasic
implicit none implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt) :: myDebug,i
integer(pInt) :: myDebug
character(len=65536) :: & character(len=pStringLen) :: &
line, & line, &
part part
character(len=pStringLen), dimension(:), allocatable :: fileContent
logical :: fileExists
write(6,'(/,a)') ' <<<+- config init -+>>>' write(6,'(/,a)') ' <<<+- config init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -120,39 +126,40 @@ subroutine config_init()
myDebug = debug_level(debug_material) myDebug = debug_level(debug_material)
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... inquire(file=trim(getSolverJobName())//'.'//material_localFileExt,exist=fileExists)
call IO_open_file(FILEUNIT,material_configFile) ! ...open material.config file if(fileExists) then
fileContent = IO_recursiveRead(trim(getSolverJobName())//'.'//material_localFileExt)
else
inquire(file='material.config',exist=fileExists)
if(.not. fileExists) call IO_error(100_pInt,ext_msg='material.config')
fileContent = IO_recursiveRead('material.config')
endif
rewind(fileUnit) do i = 1_pInt, size(fileContent)
line = '' ! to have it initialized line = trim(fileContent(i))
do while (trim(line) /= IO_EOF)
part = IO_lc(IO_getTag(line,'<','>')) part = IO_lc(IO_getTag(line,'<','>'))
select case (trim(part)) select case (trim(part))
case (trim(material_partPhase)) case (trim(material_partPhase))
call parseFile(line,phase_name,config_phase,FILEUNIT) call parseFile(phase_name,config_phase,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6)
case (trim(material_partMicrostructure)) case (trim(material_partMicrostructure))
call parseFile(line,microstructure_name,config_microstructure,FILEUNIT) call parseFile(microstructure_name,config_microstructure,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed'; flush(6)
case (trim(material_partCrystallite)) case (trim(material_partCrystallite))
call parseFile(line,crystallite_name,config_crystallite,FILEUNIT) call parseFile(crystallite_name,config_crystallite,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6)
case (trim(material_partHomogenization)) case (trim(material_partHomogenization))
call parseFile(line,homogenization_name,config_homogenization,FILEUNIT) call parseFile(homogenization_name,config_homogenization,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed'; flush(6)
case (trim(material_partTexture)) case (trim(material_partTexture))
call parseFile(line,texture_name,config_texture,FILEUNIT) call parseFile(texture_name,config_texture,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6)
case default
line = IO_read(fileUnit)
end select end select
enddo enddo
@ -173,107 +180,83 @@ end subroutine config_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief parses the material.config file !> @brief parses the material.config file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine parseFile(line,& subroutine parseFile(sectionNames,part,line, &
sectionNames,part,fileUnit) fileContent)
use prec, only: &
pStringLen
use IO, only: & use IO, only: &
IO_read, &
IO_error, & IO_error, &
IO_lc, & IO_getTag
IO_getTag, &
IO_isBlank, &
IO_stringValue, &
IO_stringPos, &
IO_EOF
implicit none implicit none
integer(pInt), intent(in) :: fileUnit character(len=64), allocatable, dimension(:), intent(out) :: sectionNames
character(len=*), dimension(:), allocatable, intent(inout) :: sectionNames
type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part
character(len=65536),intent(out) :: line character(len=pStringLen), intent(inout) :: line
character(len=pStringLen), dimension(:), intent(in) :: fileContent
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: partPosition ! position of [] tags + last line in section
integer(pInt) :: s integer(pInt) :: i, j
character(len=65536) :: devNull
character(len=64) :: tag
logical :: echo logical :: echo
echo = .false. echo = .false.
allocate(part(0))
s = 0_pInt if (allocated(part)) call IO_error(161_pInt,ext_msg=trim(line))
do while (trim(line) /= IO_EOF) ! read through sections of material part allocate(partPosition(0))
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines do i = 1_pInt, size(fileContent)
foundNextPart: if (IO_getTag(line,'<','>') /= '') then line = trim(fileContent(i))
devNull = IO_read(fileUnit, .true.) ! reset IO_read to close any recursively included files if (IO_getTag(line,'<','>') /= '') exit
exit
endif foundNextPart
nextSection: if (IO_getTag(line,'[',']') /= '') then nextSection: if (IO_getTag(line,'[',']') /= '') then
s = s + 1_pInt partPosition = [partPosition, i]
part = [part, emptyList]
tag = IO_getTag(line,'[',']')
GfortranBug86033: if (.not. allocated(sectionNames)) then
allocate(sectionNames(1),source=tag)
else GfortranBug86033
sectionNames = [sectionNames,tag]
endif GfortranBug86033
cycle cycle
endif nextSection endif nextSection
chunkPos = IO_stringPos(line) if (size(partPosition) < 1_pInt) &
tag = IO_lc(IO_stringValue(trim(line),chunkPos,1_pInt)) ! extract key echo = (trim(IO_getTag(line,'/','/')) == 'echo') .or. echo
inSection: if (s > 0_pInt) then
call part(s)%add(IO_lc(trim(line)))
else inSection
echo = (trim(tag) == '/echo/')
endif inSection
enddo enddo
if (echo) then allocate(sectionNames(size(partPosition)))
do s = 1, size(sectionNames) allocate(part(size(partPosition)))
call part(s)%show()
end do partPosition = [partPosition, i] ! needed when actually storing content
end if
do i = 1_pInt, size(partPosition) -1_pInt
sectionNames(i) = trim(adjustl(fileContent(partPosition(i))))
do j = partPosition(i) + 1_pInt, partPosition(i+1) -1_pInt
call part(i)%add(trim(adjustl(fileContent(j))))
enddo
if (echo) then
write(6,*) 'section',i, '"'//trim(sectionNames(i))//'"'
call part(i)%show()
endif
enddo
end subroutine parseFile end subroutine parseFile
!--------------------------------------------------------------------------------------------------
!> @brief deallocates the linked lists that store the content of the configuration files
!--------------------------------------------------------------------------------------------------
subroutine config_deallocate(what) subroutine config_deallocate(what)
use IO, only: & use IO, only: &
IO_error IO_error
implicit none implicit none
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
integer(pInt) :: i
select case(what) select case(trim(what))
case('material.config/phase') case('material.config/phase')
do i=1, size(config_phase)
call config_phase(i)%free
enddo
deallocate(config_phase) deallocate(config_phase)
case('material.config/microstructure') case('material.config/microstructure')
do i=1, size(config_microstructure)
call config_microstructure(i)%free
enddo
deallocate(config_microstructure) deallocate(config_microstructure)
case('material.config/crystallite') case('material.config/crystallite')
do i=1, size(config_crystallite)
call config_crystallite(i)%free
enddo
deallocate(config_crystallite) deallocate(config_crystallite)
case('material.config/homogenization') case('material.config/homogenization')
do i=1, size(config_homogenization)
call config_homogenization(i)%free
enddo
deallocate(config_homogenization) deallocate(config_homogenization)
case('material.config/texture') case('material.config/texture')
do i=1, size(config_texture)
call config_texture(i)%free
enddo
deallocate(config_texture) deallocate(config_texture)
case default case default
@ -284,11 +267,17 @@ subroutine config_deallocate(what)
end subroutine config_deallocate end subroutine config_deallocate
!##################################################################################################
! The folowing functions are part of the tPartitionedStringList object
!##################################################################################################
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief add element !> @brief add element
!> @details Adds a string together with the start/end position of chunks in this string. The new !> @details Adds a string together with the start/end position of chunks in this string. The new
!! element is added at the end of the list. Empty strings are not added. All strings are converted !! element is added at the end of the list. Empty strings are not added. All strings are converted
!! to lower case !! to lower case. The data is not stored in the new element but in the current.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine add(this,string) subroutine add(this,string)
use IO, only: & use IO, only: &
@ -299,19 +288,18 @@ subroutine add(this,string)
implicit none implicit none
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: string character(len=*), intent(in) :: string
type(tPartitionedStringList), pointer :: new, item type(tPartitionedStringList), pointer :: new, temp
if (IO_isBlank(string)) return if (IO_isBlank(string)) return
allocate(new) allocate(new)
new%string%val = IO_lc (trim(string)) temp => this
new%string%pos = IO_stringPos(trim(string)) do while (associated(temp%next))
temp => temp%next
item => this
do while (associated(item%next))
item => item%next
enddo enddo
item%next => new temp%string%val = IO_lc (trim(string))
temp%string%pos = IO_stringPos(trim(string))
temp%next => new
end subroutine add end subroutine add
@ -323,12 +311,12 @@ end subroutine add
subroutine show(this) subroutine show(this)
implicit none implicit none
class(tPartitionedStringList) :: this class(tPartitionedStringList), target, intent(in) :: this
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
item => this%next item => this
do while (associated(item)) do while (associated(item%next))
write(6,'(a)') trim(item%string%val) write(6,'(a)') ' '//trim(item%string%val)
item => item%next item => item%next
end do end do
@ -336,28 +324,55 @@ end subroutine show
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief cleans entire list !> @brief empties list and frees associated memory
!> @details list head is remains alive !> @details explicit interface to reset list. Triggers final statement (and following chain reaction)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine free(this) subroutine free(this)
implicit none implicit none
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), intent(inout) :: this
type(tPartitionedStringList), pointer :: new, item
if (.not. associated(this%next)) return if(associated(this%next)) deallocate(this%next)
item => this%next
do while (associated(item%next))
new => item
deallocate(item)
item => new%next
enddo
deallocate(item)
end subroutine free end subroutine free
!--------------------------------------------------------------------------------------------------
!> @brief empties list and frees associated memory
!> @details called when variable goes out of scope. Triggers chain reaction for list
!--------------------------------------------------------------------------------------------------
recursive subroutine finalize(this)
implicit none
type(tPartitionedStringList), intent(inout) :: this
if(associated(this%next)) deallocate(this%next)
end subroutine finalize
!--------------------------------------------------------------------------------------------------
!> @brief cleans entire array of linke lists
!> @details called when variable goes out of scope and deallocates the list at each array entry
!--------------------------------------------------------------------------------------------------
subroutine finalizeArray(this)
implicit none
integer :: i
type(tPartitionedStringList), intent(inout), dimension(:) :: this
type(tPartitionedStringList), pointer :: temp ! bug in Gfortran?
do i=1, size(this)
if (associated(this(i)%next)) then
temp => this(i)%next
!deallocate(this(i)) !internal compiler error: in gfc_build_final_call, at fortran/trans.c:975
deallocate(temp)
endif
enddo
end subroutine finalizeArray
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief reports wether a given key (string value at first position) exists in the list !> @brief reports wether a given key (string value at first position) exists in the list
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -366,14 +381,14 @@ logical function keyExists(this,key)
IO_stringValue IO_stringValue
implicit none implicit none
class(tPartitionedStringList), intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
keyExists = .false. keyExists = .false.
item => this%next item => this
do while (associated(item) .and. .not. keyExists) do while (associated(item%next) .and. .not. keyExists)
keyExists = trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key) keyExists = trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)
item => item%next item => item%next
end do end do
@ -391,14 +406,14 @@ integer(pInt) function countKeys(this,key)
implicit none implicit none
class(tPartitionedStringList), intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
countKeys = 0_pInt countKeys = 0_pInt
item => this%next item => this
do while (associated(item)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) & if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) &
countKeys = countKeys + 1_pInt countKeys = countKeys + 1_pInt
item => item%next item => item%next
@ -419,17 +434,17 @@ real(pReal) function getFloat(this,key,defaultVal)
IO_FloatValue IO_FloatValue
implicit none implicit none
class(tPartitionedStringList), intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
real(pReal), intent(in), optional :: defaultVal real(pReal), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
logical :: found logical :: found
found = present(defaultVal) found = present(defaultVal)
if (found) getFloat = defaultVal if (found) getFloat = defaultVal
item => this%next item => this
do while (associated(item)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
@ -455,17 +470,17 @@ integer(pInt) function getInt(this,key,defaultVal)
IO_IntValue IO_IntValue
implicit none implicit none
class(tPartitionedStringList), intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
integer(pInt), intent(in), optional :: defaultVal integer(pInt), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
logical :: found logical :: found
found = present(defaultVal) found = present(defaultVal)
if (found) getInt = defaultVal if (found) getInt = defaultVal
item => this%next item => this
do while (associated(item)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
@ -491,13 +506,13 @@ character(len=65536) function getString(this,key,defaultVal,raw)
IO_stringValue IO_stringValue
implicit none implicit none
class(tPartitionedStringList), intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
character(len=65536), intent(in), optional :: defaultVal character(len=65536), intent(in), optional :: defaultVal
logical, intent(in), optional :: raw logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
logical :: found, & logical :: found, &
whole whole
whole = merge(raw,.false.,present(raw)) ! whole string or white space splitting whole = merge(raw,.false.,present(raw)) ! whole string or white space splitting
found = present(defaultVal) found = present(defaultVal)
@ -506,8 +521,8 @@ character(len=65536) function getString(this,key,defaultVal,raw)
if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0_pInt,ext_msg='getString') if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0_pInt,ext_msg='getString')
endif endif
item => this%next item => this
do while (associated(item)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
@ -539,7 +554,7 @@ function getFloats(this,key,defaultVal,requiredShape)
implicit none implicit none
real(pReal), dimension(:), allocatable :: getFloats real(pReal), dimension(:), allocatable :: getFloats
class(tPartitionedStringList), intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
real(pReal), dimension(:), intent(in), optional :: defaultVal real(pReal), dimension(:), intent(in), optional :: defaultVal
integer(pInt), dimension(:), intent(in), optional :: requiredShape integer(pInt), dimension(:), intent(in), optional :: requiredShape
@ -553,8 +568,8 @@ function getFloats(this,key,defaultVal,requiredShape)
allocate(getFloats(0)) allocate(getFloats(0))
item => this%next item => this
do while (associated(item)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (.not. cumulative) getFloats = [real(pReal)::] if (.not. cumulative) getFloats = [real(pReal)::]
@ -586,7 +601,7 @@ function getInts(this,key,defaultVal,requiredShape)
implicit none implicit none
integer(pInt), dimension(:), allocatable :: getInts integer(pInt), dimension(:), allocatable :: getInts
class(tPartitionedStringList), intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
integer(pInt), dimension(:), intent(in), optional :: defaultVal, & integer(pInt), dimension(:), intent(in), optional :: defaultVal, &
requiredShape requiredShape
@ -600,8 +615,8 @@ function getInts(this,key,defaultVal,requiredShape)
allocate(getInts(0)) allocate(getInts(0))
item => this%next item => this
do while (associated(item)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (.not. cumulative) getInts = [integer(pInt)::] if (.not. cumulative) getInts = [integer(pInt)::]
@ -633,7 +648,7 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
implicit none implicit none
character(len=65536),dimension(:), allocatable :: getStrings character(len=65536),dimension(:), allocatable :: getStrings
class(tPartitionedStringList), intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
character(len=65536),dimension(:), intent(in), optional :: defaultVal character(len=65536),dimension(:), intent(in), optional :: defaultVal
integer(pInt), dimension(:), intent(in), optional :: requiredShape integer(pInt), dimension(:), intent(in), optional :: requiredShape
@ -649,8 +664,8 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
whole = merge(raw,.false.,present(raw)) whole = merge(raw,.false.,present(raw))
found = .false. found = .false.
item => this%next item => this
do while (associated(item)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings) if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings)

7
src/constitutive.f90
View File

@ -58,14 +58,15 @@ subroutine constitutive_init()
IO_write_jobIntFile, & IO_write_jobIntFile, &
IO_timeStamp IO_timeStamp
use config, only: & use config, only: &
config_deallocate config_phase
use mesh, only: & use mesh, only: &
FE_geomtype FE_geomtype
use config, only: & use config, only: &
material_Nphase, & material_Nphase, &
material_localFileExt, & material_localFileExt, &
phase_name, & phase_name, &
material_configFile material_configFile, &
config_deallocate
use material, only: & use material, only: &
material_phase, & material_phase, &
phase_plasticity, & phase_plasticity, &
@ -138,7 +139,7 @@ subroutine constitutive_init()
use kinematics_hydrogen_strain use kinematics_hydrogen_strain
implicit none implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt), parameter :: FILEUNIT = 204_pInt
integer(pInt) :: & integer(pInt) :: &
o, & !< counter in output loop o, & !< counter in output loop
ph, & !< counter in phase loop ph, & !< counter in phase loop

4
src/crystallite.f90
View File

@ -172,9 +172,9 @@ subroutine crystallite_init
IO_error IO_error
use material use material
use config, only: & use config, only: &
config_deallocate, &
config_crystallite, & config_crystallite, &
crystallite_name, & crystallite_name
config_deallocate
use constitutive, only: & use constitutive, only: &
constitutive_initialFi, & constitutive_initialFi, &
constitutive_microstructure ! derived (shortcut) quantities of given state constitutive_microstructure ! derived (shortcut) quantities of given state

2
src/debug.f90
View File

@ -102,7 +102,7 @@ subroutine debug_init
IO_EOF IO_EOF
implicit none implicit none
integer(pInt), parameter :: FILEUNIT = 300_pInt integer(pInt), parameter :: FILEUNIT = 330_pInt
integer(pInt) :: i, what integer(pInt) :: i, what
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos

114
src/homogenization.f90
View File

@ -6,9 +6,6 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module homogenization module homogenization
use prec, only: & use prec, only: &
#ifdef FEM
tOutputData, &
#endif
pInt, & pInt, &
pReal pReal
@ -22,16 +19,8 @@ module homogenization
materialpoint_P !< first P--K stress of IP materialpoint_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: &
materialpoint_dPdF !< tangent of first P--K stress at IP materialpoint_dPdF !< tangent of first P--K stress at IP
#ifdef FEM
type(tOutputData), dimension(:), allocatable, public :: &
homogOutput
type(tOutputData), dimension(:,:), allocatable, public :: &
crystalliteOutput, &
phaseOutput
#else
real(pReal), dimension(:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:), allocatable, public :: &
materialpoint_results !< results array of material point materialpoint_results !< results array of material point
#endif
integer(pInt), public, protected :: & integer(pInt), public, protected :: &
materialpoint_sizeResults, & materialpoint_sizeResults, &
homogenization_maxSizePostResults, & homogenization_maxSizePostResults, &
@ -90,20 +79,15 @@ subroutine homogenization_init
mesh_element, & mesh_element, &
FE_Nips, & FE_Nips, &
FE_geomtype FE_geomtype
#ifdef FEM
use crystallite, only: &
crystallite_sizePostResults
#else
use constitutive, only: & use constitutive, only: &
constitutive_plasticity_maxSizePostResults, & constitutive_plasticity_maxSizePostResults, &
constitutive_source_maxSizePostResults constitutive_source_maxSizePostResults
use crystallite, only: & use crystallite, only: &
crystallite_maxSizePostResults crystallite_maxSizePostResults
#endif
use config, only: & use config, only: &
config_deallocate, &
material_configFile, & material_configFile, &
material_localFileExt, & material_localFileExt, &
config_deallocate, &
config_homogenization, & config_homogenization, &
homogenization_name homogenization_name
use material use material
@ -411,33 +395,6 @@ subroutine homogenization_init
hydrogenflux_maxSizePostResults = max(hydrogenflux_maxSizePostResults ,hydrogenfluxState(p)%sizePostResults) hydrogenflux_maxSizePostResults = max(hydrogenflux_maxSizePostResults ,hydrogenfluxState(p)%sizePostResults)
enddo enddo
#ifdef FEM
allocate(homogOutput (material_Nhomogenization ))
allocate(crystalliteOutput(material_Ncrystallite, homogenization_maxNgrains))
allocate(phaseOutput (material_Nphase, homogenization_maxNgrains))
do p = 1, material_Nhomogenization
homogOutput(p)%sizeResults = homogState (p)%sizePostResults + &
thermalState (p)%sizePostResults + &
damageState (p)%sizePostResults + &
vacancyfluxState (p)%sizePostResults + &
porosityState (p)%sizePostResults + &
hydrogenfluxState(p)%sizePostResults
homogOutput(p)%sizeIpCells = count(material_homog==p)
allocate(homogOutput(p)%output(homogOutput(p)%sizeResults,homogOutput(p)%sizeIpCells))
enddo
do p = 1, material_Ncrystallite; do e = 1, homogenization_maxNgrains
crystalliteOutput(p,e)%sizeResults = crystallite_sizePostResults(p)
crystalliteOutput(p,e)%sizeIpCells = count(microstructure_crystallite(mesh_element(4,:)) == p .and. &
homogenization_Ngrains (mesh_element(3,:)) >= e)*mesh_maxNips
allocate(crystalliteOutput(p,e)%output(crystalliteOutput(p,e)%sizeResults,crystalliteOutput(p,e)%sizeIpCells))
enddo; enddo
do p = 1, material_Nphase; do e = 1, homogenization_maxNgrains
phaseOutput(p,e)%sizeResults = plasticState (p)%sizePostResults + &
sum(sourceState (p)%p(:)%sizePostResults)
phaseOutput(p,e)%sizeIpCells = count(material_phase(e,:,:) == p)
allocate(phaseOutput(p,e)%output(phaseOutput(p,e)%sizeResults,phaseOutput(p,e)%sizeIpCells))
enddo; enddo
#else
materialpoint_sizeResults = 1 & ! grain count materialpoint_sizeResults = 1 & ! grain count
+ 1 + homogenization_maxSizePostResults & ! homogSize & homogResult + 1 + homogenization_maxSizePostResults & ! homogSize & homogResult
+ thermal_maxSizePostResults & + thermal_maxSizePostResults &
@ -449,7 +406,6 @@ subroutine homogenization_init
+ 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results + 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results
+ constitutive_source_maxSizePostResults) + constitutive_source_maxSizePostResults)
allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems)) allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems))
#endif
write(6,'(/,a)') ' <<<+- homogenization init -+>>>' write(6,'(/,a)') ' <<<+- homogenization init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -473,9 +429,6 @@ subroutine homogenization_init
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
#ifndef FEM
write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_results: ', shape(materialpoint_results)
#endif
write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
endif endif
flush(6) flush(6)
@ -904,33 +857,18 @@ subroutine materialpoint_postResults
mesh_element mesh_element
use material, only: & use material, only: &
mappingHomogenization, & mappingHomogenization, &
#ifdef FEM
phaseAt, phasememberAt, &
homogenization_maxNgrains, &
material_Ncrystallite, &
material_Nphase, &
#else
homogState, & homogState, &
thermalState, & thermalState, &
damageState, & damageState, &
vacancyfluxState, & vacancyfluxState, &
porosityState, & porosityState, &
hydrogenfluxState, & hydrogenfluxState, &
#endif
plasticState, & plasticState, &
sourceState, & sourceState, &
material_phase, & material_phase, &
homogenization_Ngrains, & homogenization_Ngrains, &
microstructure_crystallite microstructure_crystallite
#ifdef FEM
use constitutive, only: &
constitutive_plasticity_maxSizePostResults, &
constitutive_source_maxSizePostResults
#endif
use crystallite, only: & use crystallite, only: &
#ifdef FEM
crystallite_maxSizePostResults, &
#endif
crystallite_sizePostResults, & crystallite_sizePostResults, &
crystallite_postResults crystallite_postResults
@ -943,55 +881,6 @@ subroutine materialpoint_postResults
g, & !< grain number g, & !< grain number
i, & !< integration point number i, & !< integration point number
e !< element number e !< element number
#ifdef FEM
integer(pInt) :: &
myHomog, &
myPhase, &
crystalliteCtr(material_Ncrystallite, homogenization_maxNgrains), &
phaseCtr (material_Nphase, homogenization_maxNgrains)
real(pReal), dimension(1+crystallite_maxSizePostResults + &
1+constitutive_plasticity_maxSizePostResults + &
constitutive_source_maxSizePostResults) :: &
crystalliteResults
crystalliteCtr = 0_pInt; phaseCtr = 0_pInt
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
myCrystallite = microstructure_crystallite(mesh_element(4,e))
IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
myHomog = mappingHomogenization(2,i,e)
thePos = mappingHomogenization(1,i,e)
homogOutput(myHomog)%output(1: &
homogOutput(myHomog)%sizeResults, &
thePos) = homogenization_postResults(i,e)
grainLooping :do g = 1,myNgrains
myPhase = phaseAt(g,i,e)
crystalliteResults(1:1+crystallite_sizePostResults(myCrystallite) + &
1+plasticState(myPhase)%sizePostResults + &
sum(sourceState(myPhase)%p(:)%sizePostResults)) = crystallite_postResults(g,i,e)
if (microstructure_crystallite(mesh_element(4,e)) == myCrystallite .and. &
homogenization_Ngrains (mesh_element(3,e)) >= g) then
crystalliteCtr(myCrystallite,g) = crystalliteCtr(myCrystallite,g) + 1_pInt
crystalliteOutput(myCrystallite,g)% &
output(1:crystalliteOutput(myCrystallite,g)%sizeResults,crystalliteCtr(myCrystallite,g)) = &
crystalliteResults(2:1+crystalliteOutput(myCrystallite,g)%sizeResults)
endif
if (material_phase(g,i,e) == myPhase) then
phaseCtr(myPhase,g) = phaseCtr(myPhase,g) + 1_pInt
phaseOutput(myPhase,g)% &
output(1:phaseOutput(myPhase,g)%sizeResults,phaseCtr(myPhase,g)) = &
crystalliteResults(3 + crystalliteOutput(myCrystallite,g)%sizeResults: &
1 + crystalliteOutput(myCrystallite,g)%sizeResults + &
1 + plasticState (myphase)%sizePostResults + &
sum(sourceState(myphase)%p(:)%sizePostResults))
endif
enddo grainLooping
enddo IpLooping
enddo elementLooping
#else
!$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,thePos,theSize) !$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,thePos,theSize)
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -1027,7 +916,6 @@ subroutine materialpoint_postResults
enddo IpLooping enddo IpLooping
enddo elementLooping enddo elementLooping
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
#endif
end subroutine materialpoint_postResults end subroutine materialpoint_postResults

48
src/lattice.f90
View File

@ -111,6 +111,9 @@ module lattice
-1,-1, 0, -1, 1,-1 & ! D6 -1,-1, 0, -1, 1,-1 & ! D6
],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Nslip]) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli ],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Nslip]) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli
character(len=*), dimension(1), parameter, public :: LATTICE_FCC_SLIPFAMILY_NAME = &
['<0 1 -1>{1 1 1}']
real(pReal), dimension(3+3,LATTICE_fcc_Ntwin), parameter, private :: & real(pReal), dimension(3+3,LATTICE_fcc_Ntwin), parameter, private :: &
LATTICE_fcc_systemTwin = reshape(real( [& LATTICE_fcc_systemTwin = reshape(real( [&
-2, 1, 1, 1, 1, 1, & -2, 1, 1, 1, 1, 1, &
@ -127,6 +130,9 @@ module lattice
-1, 1, 2, -1, 1,-1 & -1, 1, 2, -1, 1,-1 &
],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ntwin]) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli ],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ntwin]) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli
character(len=*), dimension(1), parameter, public :: LATTICE_FCC_TWINFAMILY_NAME = &
['<-2 1 1>{1 1 1}']
real(pReal), dimension(3+3,LATTICE_fcc_Ntrans), parameter, private :: & real(pReal), dimension(3+3,LATTICE_fcc_Ntrans), parameter, private :: &
LATTICE_fccTohex_systemTrans = reshape(real( [& LATTICE_fccTohex_systemTrans = reshape(real( [&
-2, 1, 1, 1, 1, 1, & -2, 1, 1, 1, 1, 1, &
@ -433,6 +439,10 @@ module lattice
! 1,-1, 1, 3, 2,-1 & ! 1,-1, 1, 3, 2,-1 &
],pReal),[ 3_pInt + 3_pInt ,LATTICE_bcc_Nslip]) ],pReal),[ 3_pInt + 3_pInt ,LATTICE_bcc_Nslip])
character(len=*), dimension(2), parameter, public :: LATTICE_BCC_SLIPFAMILY_NAME = &
['<1 -1 1>{0 1 1}', &
'<1 -1 1>{2 1 1}']
real(pReal), dimension(3+3,LATTICE_bcc_Ntwin), parameter, private :: & real(pReal), dimension(3+3,LATTICE_bcc_Ntwin), parameter, private :: &
LATTICE_bcc_systemTwin = reshape(real([& LATTICE_bcc_systemTwin = reshape(real([&
! Twin system <111>{112} ! Twin system <111>{112}
@ -450,6 +460,9 @@ module lattice
1, 1, 1, 1, 1,-2 & 1, 1, 1, 1, 1,-2 &
],pReal),[ 3_pInt + 3_pInt,LATTICE_bcc_Ntwin]) ],pReal),[ 3_pInt + 3_pInt,LATTICE_bcc_Ntwin])
character(len=*), dimension(1), parameter, public :: LATTICE_BCC_TWINFAMILY_NAME = &
['<1 1 1>{2 1 1}']
real(pReal), dimension(LATTICE_bcc_Ntwin), parameter, private :: & real(pReal), dimension(LATTICE_bcc_Ntwin), parameter, private :: &
LATTICE_bcc_shearTwin = 0.5_pReal*sqrt(2.0_pReal) LATTICE_bcc_shearTwin = 0.5_pReal*sqrt(2.0_pReal)
@ -618,6 +631,14 @@ module lattice
1, 1, -2, 3, -1, -1, 2, 2 & 1, 1, -2, 3, -1, -1, 2, 2 &
],pReal),[ 4_pInt + 4_pInt,LATTICE_hex_Nslip]) !< slip systems for hex sorted by A. Alankar & P. Eisenlohr ],pReal),[ 4_pInt + 4_pInt,LATTICE_hex_Nslip]) !< slip systems for hex sorted by A. Alankar & P. Eisenlohr
character(len=*), dimension(6), parameter, public :: LATTICE_HEX_SLIPFAMILY_NAME = &
['<1 1 . 1>{0 0 . 1} ', &
'<1 1 . 1>{1 0 . 0} ', &
'<1 0 . 0>{1 1 . 0} ', &
'<1 1 . 0>{-1 1 . 1} ', &
'<1 1 . 3>{-1 0 . 1} ', &
'<1 1 . 3>{-1 -1 . 2}']
real(pReal), dimension(4+4,LATTICE_hex_Ntwin), parameter, private :: & real(pReal), dimension(4+4,LATTICE_hex_Ntwin), parameter, private :: &
LATTICE_hex_systemTwin = reshape(real([& LATTICE_hex_systemTwin = reshape(real([&
! Compression or Tension =f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981) ! Compression or Tension =f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981)
@ -650,6 +671,12 @@ module lattice
1, 1, -2, -3, 1, 1, -2, 2 & 1, 1, -2, -3, 1, 1, -2, 2 &
],pReal),[ 4_pInt + 4_pInt ,LATTICE_hex_Ntwin]) !< twin systems for hex, order follows Prof. Tom Bieler's scheme; but numbering in data was restarted from 1 ],pReal),[ 4_pInt + 4_pInt ,LATTICE_hex_Ntwin]) !< twin systems for hex, order follows Prof. Tom Bieler's scheme; but numbering in data was restarted from 1
character(len=*), dimension(4), parameter, public :: LATTICE_HEX_TWINFAMILY_NAME = &
['<-1 0 . 1>{1 0 . 2} ', &
'<1 1 . 6>{-1 -1 . 1}', &
'<1 0 . -2>{1 0 . 1} ', &
'<1 1 . -3>{1 1 . 2} ']
integer(pInt), dimension(LATTICE_hex_Ntwin), parameter, private :: & integer(pInt), dimension(LATTICE_hex_Ntwin), parameter, private :: &
LATTICE_hex_shearTwin = reshape(int( [& ! indicator to formula further below LATTICE_hex_shearTwin = reshape(int( [& ! indicator to formula further below
1, & ! <-10.1>{10.2} 1, & ! <-10.1>{10.2}
@ -926,6 +953,21 @@ module lattice
1, 1, 1, 1,-2, 1 & 1, 1, 1, 1,-2, 1 &
],pReal),[ 3_pInt + 3_pInt,LATTICE_bct_Nslip]) !< slip systems for bct sorted by Bieler ],pReal),[ 3_pInt + 3_pInt,LATTICE_bct_Nslip]) !< slip systems for bct sorted by Bieler
character(len=*), dimension(13), parameter, public :: LATTICE_BCT_SLIPFAMILY_NAME = &
['{1 0 0)<0 0 1] ', &
'{1 1 0)<0 0 1] ', &
'{1 0 0)<0 1 0] ', &
'{1 1 0)<1 -1 1]', &
'{1 1 0)<1 -1 0]', &
'{1 0 0)<0 1 1] ', &
'{0 0 1)<0 1 0] ', &
'{0 0 1)<1 1 0] ', &
'{0 1 1)<0 1 -1]', &
'{0 1 1)<1 -1 1]', &
'{0 1 1)<1 0 0] ', &
'{2 1 1)<0 1 -1]', &
'{2 1 1)<-1 1 1]']
integer(pInt), dimension(LATTICE_bct_Nslip,LATTICE_bct_Nslip), parameter, public :: & integer(pInt), dimension(LATTICE_bct_Nslip,LATTICE_bct_Nslip), parameter, public :: &
LATTICE_bct_interactionSlipSlip = reshape(int( [& LATTICE_bct_interactionSlipSlip = reshape(int( [&
1, 2, 3, 3, 7, 7, 13, 13, 13, 13, 21, 21, 31, 31, 31, 31, 43, 43, 57, 57, 73, 73, 73, 73, 91, 91, 91, 91, 91, 91, 91, 91, 111, 111, 111, 111, 133,133,133,133,133,133,133,133, 157,157,157,157,157,157,157,157, & 1, 2, 3, 3, 7, 7, 13, 13, 13, 13, 21, 21, 31, 31, 31, 31, 43, 43, 57, 57, 73, 73, 73, 73, 91, 91, 91, 91, 91, 91, 91, 91, 111, 111, 111, 111, 133,133,133,133,133,133,133,133, 157,157,157,157,157,157,157,157, &
@ -1278,7 +1320,7 @@ subroutine lattice_init
integer(pInt) :: Nphases integer(pInt) :: Nphases
character(len=65536) :: & character(len=65536) :: &
tag = '' tag = ''
integer(pInt) :: section = 0_pInt,i,p integer(pInt) :: i,p
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
temp, & temp, &
CoverA, & !< c/a ratio for low symmetry type lattice CoverA, & !< c/a ratio for low symmetry type lattice
@ -1388,9 +1430,9 @@ subroutine lattice_init
tag = config_phase(p)%getString('trans_lattice_structure',defaultVal=tag) tag = config_phase(p)%getString('trans_lattice_structure',defaultVal=tag)
select case(trim(tag)) select case(trim(tag))
case('bcc') case('bcc')
trans_lattice_structure(section) = LATTICE_bcc_ID trans_lattice_structure(p) = LATTICE_bcc_ID
case('hex','hexagonal') case('hex','hexagonal')
trans_lattice_structure(section) = LATTICE_hex_ID trans_lattice_structure(p) = LATTICE_hex_ID
end select end select
lattice_C66(1,1,p) = config_phase(p)%getFloat('c11',defaultVal=0.0_pReal) lattice_C66(1,1,p) = config_phase(p)%getFloat('c11',defaultVal=0.0_pReal)

16
src/material.f90
View File

@ -16,7 +16,7 @@ module material
tSourceState, & tSourceState, &
tHomogMapping, & tHomogMapping, &
tPhaseMapping, & tPhaseMapping, &
group_scalar, & group_float, &
group_int group_int
implicit none implicit none
@ -268,7 +268,7 @@ module material
porosityMapping, & !< mapping for porosity state/fields porosityMapping, & !< mapping for porosity state/fields
hydrogenfluxMapping !< mapping for hydrogen conc state/fields hydrogenfluxMapping !< mapping for hydrogen conc state/fields
type(group_scalar), allocatable, dimension(:), public :: & type(group_float), allocatable, dimension(:), public :: &
temperature, & !< temperature field temperature, & !< temperature field
damage, & !< damage field damage, & !< damage field
vacancyConc, & !< vacancy conc field vacancyConc, & !< vacancy conc field
@ -360,8 +360,7 @@ subroutine material_init()
homogenization_name, & homogenization_name, &
microstructure_name, & microstructure_name, &
phase_name, & phase_name, &
texture_name, & texture_name
config_deallocate
use mesh, only: & use mesh, only: &
mesh_maxNips, & mesh_maxNips, &
mesh_NcpElems, & mesh_NcpElems, &
@ -370,7 +369,7 @@ subroutine material_init()
FE_geomtype FE_geomtype
implicit none implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt), parameter :: FILEUNIT = 210_pInt
integer(pInt) :: m,c,h, myDebug, myPhase, myHomog integer(pInt) :: m,c,h, myDebug, myPhase, myHomog
integer(pInt) :: & integer(pInt) :: &
g, & !< grain number g, & !< grain number
@ -469,7 +468,6 @@ subroutine material_init()
endif debugOut endif debugOut
call material_populateGrains call material_populateGrains
call config_deallocate('material.config/microstructure')
allocate(phaseAt ( homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt) allocate(phaseAt ( homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
allocate(phasememberAt ( homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt) allocate(phasememberAt ( homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
@ -921,8 +919,8 @@ subroutine material_parseTexture
IO_floatValue, & IO_floatValue, &
IO_stringValue IO_stringValue
use config, only: & use config, only: &
config_texture, & config_deallocate, &
config_deallocate config_texture
use math, only: & use math, only: &
inRad, & inRad, &
math_sampleRandomOri, & math_sampleRandomOri, &
@ -1093,6 +1091,7 @@ subroutine material_populateGrains
use config, only: & use config, only: &
config_homogenization, & config_homogenization, &
config_microstructure, & config_microstructure, &
config_deallocate, &
homogenization_name, & homogenization_name, &
microstructure_name microstructure_name
use IO, only: & use IO, only: &
@ -1429,6 +1428,7 @@ subroutine material_populateGrains
deallocate(texture_transformation) deallocate(texture_transformation)
deallocate(Nelems) deallocate(Nelems)
deallocate(elemsOfHomogMicro) deallocate(elemsOfHomogMicro)
call config_deallocate('material.config/microstructure')
end subroutine material_populateGrains end subroutine material_populateGrains

10
src/math.f90
View File

@ -12,7 +12,7 @@ module math
implicit none implicit none
private private
real(pReal), parameter, public :: PI = 3.141592653589793_pReal !< ratio of a circle's circumference to its diameter real(pReal), parameter, public :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter
real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree
real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian
complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)*(PI,0.0_pReal) !< Re(0.0), Im(2xPi) complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)*(PI,0.0_pReal) !< Re(0.0), Im(2xPi)
@ -36,13 +36,13 @@ module math
real(pReal), dimension(6), parameter, private :: & real(pReal), dimension(6), parameter, private :: &
nrmMandel = [& nrmMandel = [&
1.0_pReal, 1.0_pReal, 1.0_pReal,& 1.0_pReal, 1.0_pReal, 1.0_pReal, &
1.414213562373095_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal ] !< weighting for Mandel notation (forward) sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< weighting for Mandel notation (forward)
real(pReal), dimension(6), parameter , public :: & real(pReal), dimension(6), parameter , public :: &
invnrmMandel = [& invnrmMandel = [&
1.0_pReal, 1.0_pReal, 1.0_pReal,& 1.0_pReal, 1.0_pReal, 1.0_pReal, &
0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal ] !< weighting for Mandel notation (backward) 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward)
integer(pInt), dimension (2,6), parameter, private :: & integer(pInt), dimension (2,6), parameter, private :: &
mapVoigt = reshape([& mapVoigt = reshape([&

33
src/mesh.f90
View File

@ -95,9 +95,11 @@ module mesh
integer(pInt), dimension(:,:), allocatable, private :: & integer(pInt), dimension(:,:), allocatable, private :: &
mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID
#if defined(Marc4DAMASK) || defined(Abaqus)
integer(pInt), dimension(:,:), allocatable, target, private :: & integer(pInt), dimension(:,:), allocatable, target, private :: &
mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid]
mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid] mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid]
#endif
integer(pInt),dimension(:,:,:), allocatable, private :: & integer(pInt),dimension(:,:,:), allocatable, private :: &
mesh_cell !< cell connectivity for each element,ip/cell mesh_cell !< cell connectivity for each element,ip/cell
@ -402,7 +404,9 @@ module mesh
public :: & public :: &
mesh_init, & mesh_init, &
#if defined(Marc4DAMASK) || defined(Abaqus)
mesh_FEasCP, & mesh_FEasCP, &
#endif
mesh_build_cellnodes, & mesh_build_cellnodes, &
mesh_build_ipVolumes, & mesh_build_ipVolumes, &
mesh_build_ipCoordinates, & mesh_build_ipCoordinates, &
@ -420,7 +424,6 @@ module mesh
#ifdef Spectral #ifdef Spectral
mesh_spectral_getHomogenization, & mesh_spectral_getHomogenization, &
mesh_spectral_count, & mesh_spectral_count, &
mesh_spectral_mapNodesAndElems, &
mesh_spectral_count_cpSizes, & mesh_spectral_count_cpSizes, &
mesh_spectral_build_nodes, & mesh_spectral_build_nodes, &
mesh_spectral_build_elements, & mesh_spectral_build_elements, &
@ -552,8 +555,6 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6) if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6)
call mesh_spectral_count() call mesh_spectral_count()
if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6) if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6)
call mesh_spectral_mapNodesAndElems
if (myDebug) write(6,'(a)') ' Mapped nodes and elements'; flush(6)
call mesh_spectral_count_cpSizes call mesh_spectral_count_cpSizes
if (myDebug) write(6,'(a)') ' Built CP statistics'; flush(6) if (myDebug) write(6,'(a)') ' Built CP statistics'; flush(6)
call mesh_spectral_build_nodes() call mesh_spectral_build_nodes()
@ -659,12 +660,16 @@ subroutine mesh_init(ip,el)
allocate(calcMode(mesh_maxNips,mesh_NcpElems)) allocate(calcMode(mesh_maxNips,mesh_NcpElems))
calcMode = .false. ! pretend to have collected what first call is asking (F = I) calcMode = .false. ! pretend to have collected what first call is asking (F = I)
#if defined(Marc4DAMASK) || defined(Abaqus)
calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc"
#else
calcMode(ip,el) = .true. ! first ip,el needs to be already pingponged to "calc"
#endif
end subroutine mesh_init end subroutine mesh_init
#if defined(Marc4DAMASK) || defined(Abaqus)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Gives the FE to CP ID mapping by binary search through lookup array !> @brief Gives the FE to CP ID mapping by binary search through lookup array
!! valid questions (what) are 'elem', 'node' !! valid questions (what) are 'elem', 'node'
@ -713,7 +718,7 @@ integer(pInt) function mesh_FEasCP(what,myID)
enddo binarySearch enddo binarySearch
end function mesh_FEasCP end function mesh_FEasCP
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Split CP elements into cells. !> @brief Split CP elements into cells.
@ -1188,24 +1193,6 @@ subroutine mesh_spectral_count()
end subroutine mesh_spectral_count end subroutine mesh_spectral_count
!--------------------------------------------------------------------------------------------------
!> @brief fake map node from FE ID to internal (consecutive) representation for node and element
!! Allocates global array 'mesh_mapFEtoCPnode' and 'mesh_mapFEtoCPelem'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_mapNodesAndElems
use math, only: &
math_range
implicit none
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes), source = 0_pInt)
allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems), source = 0_pInt)
mesh_mapFEtoCPnode = spread(math_range(mesh_Nnodes),1,2)
mesh_mapFEtoCPelem = spread(math_range(mesh_NcpElems),1,2)
end subroutine mesh_spectral_mapNodesAndElems
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements. !> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements.
!! Sets global values 'mesh_maxNnodes', 'mesh_maxNips', 'mesh_maxNipNeighbors', !! Sets global values 'mesh_maxNnodes', 'mesh_maxNips', 'mesh_maxNipNeighbors',

356
src/meshFEM.f90 Normal file
View File

@ -0,0 +1,356 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Driver controlling inner and outer load case looping of the FEM solver
!> @details doing cutbacking, forwarding in case of restart, reporting statistics, writing
!> results
!--------------------------------------------------------------------------------------------------
module mesh
#include <petsc/finclude/petscis.h>
#include <petsc/finclude/petscdmda.h>
use prec, only: pReal, pInt
use PETScdmda
use PETScis
implicit none
private
integer(pInt), public, protected :: &
mesh_Nboundaries, &
mesh_NcpElems, & !< total number of CP elements in mesh
mesh_NcpElemsGlobal, &
mesh_Nnodes, & !< total number of nodes in mesh
mesh_maxNnodes, & !< max number of nodes in any CP element
mesh_maxNips, & !< max number of IPs in any CP element
mesh_maxNipNeighbors, &
mesh_Nelems !< total number of elements in mesh
real(pReal), public, protected :: charLength
integer(pInt), dimension(:,:), allocatable, public, protected :: &
mesh_element !< FEid, type(internal representation), material, texture, node indices as CP IDs
real(pReal), dimension(:,:), allocatable, public :: &
mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:), allocatable, public, protected :: &
mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!)
real(pReal), dimension(:,:,:), allocatable, public :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
mesh_ipArea !< area of interface to neighboring IP (initially!)
real(pReal),dimension(:,:,:,:), allocatable, public, protected :: &
mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!)
integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: &
mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
logical, dimension(3), public, protected :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes)
DM, public :: geomMesh
integer(pInt), dimension(:), allocatable, public, protected :: &
mesh_boundaries
integer(pInt), parameter, public :: &
FE_Nelemtypes = 1_pInt, &
FE_Ngeomtypes = 1_pInt, &
FE_Ncelltypes = 1_pInt, &
FE_maxNnodes = 1_pInt, &
FE_maxNips = 14_pInt
integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_geomtype = & !< geometry type of particular element type
int([1],pInt)
integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_celltype = & !< cell type that is used by each geometry type
int([1],pInt)
integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = & !< number of nodes that constitute a specific type of element
int([0],pInt)
integer(pInt), dimension(FE_Ngeomtypes), public :: FE_Nips = & !< number of IPs in a specific type of element
int([0],pInt)
integer(pInt), dimension(FE_Ncelltypes), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type
int([6],pInt)
public :: &
mesh_init, &
mesh_FEM_build_ipVolumes, &
mesh_FEM_build_ipCoordinates, &
mesh_cellCenterCoordinates
external :: &
DMPlexCreateFromFile, &
DMPlexDistribute, &
DMPlexCopyCoordinates, &
DMGetStratumSize, &
DMPlexGetHeightStratum, &
DMPlexGetLabelValue, &
DMPlexSetLabelValue
contains
!--------------------------------------------------------------------------------------------------
!> @brief initializes the mesh by calling all necessary private routines the mesh module
!! Order and routines strongly depend on type of solver
!--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,el)
use DAMASK_interface
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use IO, only: &
IO_timeStamp, &
IO_error, &
IO_open_file, &
IO_stringPos, &
IO_intValue, &
IO_EOF, &
IO_read, &
IO_isBlank
use debug, only: &
debug_e, &
debug_i
use numerics, only: &
usePingPong, &
integrationOrder, &
worldrank, &
worldsize
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP, &
calcMode
use FEM_Zoo, only: &
FEM_Zoo_nQuadrature, &
FEM_Zoo_QuadraturePoints
implicit none
integer(pInt), parameter :: FILEUNIT = 222_pInt
integer(pInt), intent(in) :: el, ip
integer(pInt) :: j
integer(pInt), allocatable, dimension(:) :: chunkPos
integer :: dimPlex
character(len=512) :: &
line
logical :: flag
PetscSF :: sf
DM :: globalMesh
PetscInt :: face, nFaceSets
PetscInt, pointer :: pFaceSets(:)
IS :: faceSetIS
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
CHKERRQ(ierr)
call DMGetDimension(globalMesh,dimPlex,ierr)
CHKERRQ(ierr)
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
CHKERRQ(ierr)
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
CHKERRQ(ierr)
call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
allocate(mesh_boundaries(mesh_Nboundaries), source = 0_pInt)
call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr)
CHKERRQ(ierr)
call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,ierr)
CHKERRQ(ierr)
if (nFaceSets > 0) call ISGetIndicesF90(faceSetIS,pFaceSets,ierr)
do face = 1, nFaceSets
mesh_boundaries(face) = pFaceSets(face)
enddo
if (nFaceSets > 0) call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr)
call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
if (worldrank == 0) then
j = 0
flag = .false.
call IO_open_file(FILEUNIT,trim(geometryFile))
do
read(FILEUNIT,'(a512)') line
if (trim(line) == IO_EOF) exit ! skip empty lines
if (trim(line) == '$Elements') then
read(FILEUNIT,'(a512)') line
read(FILEUNIT,'(a512)') line
flag = .true.
endif
if (trim(line) == '$EndElements') exit
if (flag) then
chunkPos = IO_stringPos(line)
if (chunkPos(1) == 3+IO_intValue(line,chunkPos,3)+dimPlex+1) then
call DMSetLabelValue(globalMesh,'material',j,IO_intValue(line,chunkPos,4),ierr)
CHKERRQ(ierr)
j = j + 1
endif ! count all identifiers to allocate memory and do sanity check
endif
enddo
close (FILEUNIT)
endif
if (worldsize > 1) then
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr)
else
call DMClone(globalMesh,geomMesh,ierr)
CHKERRQ(ierr)
endif
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_Nelems,ierr)
CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
CHKERRQ(ierr)
mesh_NcpElems = mesh_Nelems
FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
mesh_maxNnodes = FE_Nnodes(1_pInt)
mesh_maxNips = FE_Nips(1_pInt)
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
call mesh_FEM_build_ipVolumes(dimPlex)
allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)); mesh_element = 0_pInt
do j = 1, mesh_NcpElems
mesh_element( 1,j) = j
mesh_element( 2,j) = 1_pInt ! elem type
mesh_element( 3,j) = 1_pInt ! homogenization
call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
CHKERRQ(ierr)
end do
if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
call IO_error(600_pInt) ! ping-pong must be disabled when having non-DAMASK elements
if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
call IO_error(602_pInt,ext_msg='element') ! selected element does not exist
if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2_pInt,debug_e)))) &
call IO_error(602_pInt,ext_msg='IP') ! selected element does not have requested IP
FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt ! parallel loop bounds set to comprise from first IP...
forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element
if (allocated(calcMode)) deallocate(calcMode)
allocate(calcMode(mesh_maxNips,mesh_NcpElems))
calcMode = .false. ! pretend to have collected what first call is asking (F = I)
calcMode(ip,el) = .true. ! first ip,el needs to be already pingponged to "calc"
end subroutine mesh_init
!--------------------------------------------------------------------------------------------------
!> @brief Calculates cell center coordinates.
!--------------------------------------------------------------------------------------------------
pure function mesh_cellCenterCoordinates(ip,el)
implicit none
integer(pInt), intent(in) :: el, & !< element number
ip !< integration point number
real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell
end function mesh_cellCenterCoordinates
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
!> @details The IP volume is calculated differently depending on the cell type.
!> 2D cells assume an element depth of one in order to calculate the volume.
!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal
!> shape with a cell face as basis and the central ip at the tip. This subvolume is
!> calculated as an average of four tetrahedals with three corners on the cell face
!> and one corner at the central ip.
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipVolumes(dimPlex)
use math, only: &
math_I3, &
math_det33
implicit none
PetscInt :: dimPlex
PetscReal :: vol
PetscReal, target :: cent(dimPlex), norm(dimPlex)
PetscReal, pointer :: pCent(:), pNorm(:)
PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
endif
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
pCent => cent
pNorm => norm
do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr)
mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
enddo
end subroutine mesh_FEM_build_ipVolumes
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates'
! Called by all solvers in mesh_init in order to initialize the ip coordinates.
! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus,
! so no need to use this subroutine anymore; Marc however only provides nodal displacements,
! so in this case the ip coordinates are always calculated on the basis of this subroutine.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES,
! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME.
! HAS TO BE CHANGED IN A LATER VERSION.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
implicit none
PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
qOffset = 0
do qPt = 1, mesh_maxNips
do dirI = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
do dirJ = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0)
enddo
enddo
qOffset = qOffset + dimPlex
enddo
enddo
end subroutine mesh_FEM_build_ipCoordinates
end module mesh

3
src/plastic_isotropic.f90
View File

@ -109,11 +109,9 @@ use IO
type(tParameters), pointer :: prm type(tParameters), pointer :: prm
integer(pInt) :: & integer(pInt) :: &
o, &
phase, & phase, &
instance, & instance, &
maxNinstance, & maxNinstance, &
mySize, &
sizeDotState, & sizeDotState, &
sizeState, & sizeState, &
sizeDeltaState sizeDeltaState
@ -136,7 +134,6 @@ use IO
plastic_isotropic_output = '' plastic_isotropic_output = ''
allocate(plastic_isotropic_Noutput(maxNinstance), source=0_pInt) allocate(plastic_isotropic_Noutput(maxNinstance), source=0_pInt)
! inernal variable
allocate(param(maxNinstance)) ! one container of parameters per instance allocate(param(maxNinstance)) ! one container of parameters per instance
allocate(state(maxNinstance)) ! internal state aliases allocate(state(maxNinstance)) ! internal state aliases
allocate(dotState(maxNinstance)) allocate(dotState(maxNinstance))

16
src/plastic_phenopowerlaw.f90
View File

@ -241,29 +241,29 @@ subroutine plastic_phenopowerlaw_init
select case(outputs(i)) select case(outputs(i))
case ('resistance_slip') case ('resistance_slip')
outputID = resistance_slip_ID outputID = resistance_slip_ID
outputSize = sum(prm%Nslip) outputSize = prm%totalNslip
case ('accumulatedshear_slip') case ('accumulatedshear_slip')
outputID = accumulatedshear_slip_ID outputID = accumulatedshear_slip_ID
outputSize = sum(prm%Nslip) outputSize = prm%totalNslip
case ('shearrate_slip') case ('shearrate_slip')
outputID = shearrate_slip_ID outputID = shearrate_slip_ID
outputSize = sum(prm%Nslip) outputSize = prm%totalNslip
case ('resolvedstress_slip') case ('resolvedstress_slip')
outputID = resolvedstress_slip_ID outputID = resolvedstress_slip_ID
outputSize = sum(prm%Nslip) outputSize = prm%totalNslip
case ('resistance_twin') case ('resistance_twin')
outputID = resistance_twin_ID outputID = resistance_twin_ID
outputSize = sum(prm%Ntwin) outputSize = prm%totalNtwin
case ('accumulatedshear_twin') case ('accumulatedshear_twin')
outputID = accumulatedshear_twin_ID outputID = accumulatedshear_twin_ID
outputSize = sum(prm%Ntwin) outputSize = prm%totalNtwin
case ('shearrate_twin') case ('shearrate_twin')
outputID = shearrate_twin_ID outputID = shearrate_twin_ID
outputSize = sum(prm%Ntwin) outputSize = prm%totalNtwin
case ('resolvedstress_twin') case ('resolvedstress_twin')
outputID = resolvedstress_twin_ID outputID = resolvedstress_twin_ID
outputSize = sum(prm%Ntwin) outputSize = prm%totalNtwin
case ('totalvolfrac_twin') case ('totalvolfrac_twin')
outputID = totalvolfrac_twin_ID outputID = totalvolfrac_twin_ID

6
src/prec.f90
View File

@ -7,6 +7,7 @@
!> @brief setting precision for real and int type !> @brief setting precision for real and int type
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module prec module prec
! ToDo: use, intrinsic :: iso_fortran_env, only : I8 => int64, WP => real64
implicit none implicit none
private private
#if (FLOAT==8) #if (FLOAT==8)
@ -23,14 +24,15 @@ module prec
NO SUITABLE PRECISION FOR INTEGER SELECTED, STOPPING COMPILATION NO SUITABLE PRECISION FOR INTEGER SELECTED, STOPPING COMPILATION
#endif #endif
integer, parameter, public :: pStringLen = 256 !< default string lenth
integer, parameter, public :: pLongInt = 8 !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12) integer, parameter, public :: pLongInt = 8 !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12)
real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation)
integer(pInt), allocatable, dimension(:) :: realloc_lhs_test integer(pInt), allocatable, dimension(:) :: realloc_lhs_test
type, public :: group_scalar !< variable length datatype used for storage of state type, public :: group_float !< variable length datatype used for storage of state
real(pReal), dimension(:), pointer :: p real(pReal), dimension(:), pointer :: p
end type group_scalar end type group_float
type, public :: group_int type, public :: group_int
integer(pInt), dimension(:), pointer :: p integer(pInt), dimension(:), pointer :: p

56
src/system_routines.f90
View File

@ -78,28 +78,30 @@ end function isDirectory
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief gets the current working directory !> @brief gets the current working directory
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical function getCWD(str) character(len=1024) function getCWD()
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_INT, & C_INT, &
C_CHAR, & C_CHAR, &
C_NULL_CHAR C_NULL_CHAR
implicit none implicit none
character(len=*), intent(out) :: str character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array
character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string is an array
integer(C_INT) :: stat integer(C_INT) :: stat
integer :: i integer :: i
str = repeat('',len(str)) call getCurrentWorkDir_C(charArray,stat)
call getCurrentWorkDir_C(strFixedLength,stat) if (stat /= 0_C_INT) then
do i=1,1024 ! copy array components until Null string is found getCWD = 'Error occured when getting currend working directory'
if (strFixedLength(i) /= C_NULL_CHAR) then else
str(i:i)=strFixedLength(i) getCWD = repeat('',len(getCWD))
else arrayToString: do i=1,len(getCWD)
exit if (charArray(i) /= C_NULL_CHAR) then
endif getCWD(i:i)=charArray(i)
enddo else
getCWD=merge(.True.,.False.,stat /= 0_C_INT) exit
endif
enddo arrayToString
endif
end function getCWD end function getCWD
@ -107,28 +109,30 @@ end function getCWD
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief gets the current host name !> @brief gets the current host name
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical function getHostName(str) character(len=1024) function getHostName()
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_INT, & C_INT, &
C_CHAR, & C_CHAR, &
C_NULL_CHAR C_NULL_CHAR
implicit none implicit none
character(len=*), intent(out) :: str character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array
character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string is an array
integer(C_INT) :: stat integer(C_INT) :: stat
integer :: i integer :: i
str = repeat('',len(str)) call getHostName_C(charArray,stat)
call getHostName_C(strFixedLength,stat) if (stat /= 0_C_INT) then
do i=1,1024 ! copy array components until Null string is found getHostName = 'Error occured when getting host name'
if (strFixedLength(i) /= C_NULL_CHAR) then else
str(i:i)=strFixedLength(i) getHostName = repeat('',len(getHostName))
else arrayToString: do i=1,len(getHostName)
exit if (charArray(i) /= C_NULL_CHAR) then
endif getHostName(i:i)=charArray(i)
enddo else
getHostName=merge(.True.,.False.,stat /= 0_C_INT) exit
endif
enddo arrayToString
endif
end function getHostName end function getHostName

4
src/vacancyflux_cahnhilliard.f90
View File

@ -7,7 +7,7 @@ module vacancyflux_cahnhilliard
use prec, only: & use prec, only: &
pReal, & pReal, &
pInt, & pInt, &
group_scalar group_float
implicit none implicit none
private private
@ -26,7 +26,7 @@ module vacancyflux_cahnhilliard
real(pReal), dimension(:), allocatable, private :: & real(pReal), dimension(:), allocatable, private :: &
vacancyflux_cahnhilliard_flucAmplitude vacancyflux_cahnhilliard_flucAmplitude
type(group_scalar), dimension(:), allocatable, private :: & type(group_float), dimension(:), allocatable, private :: &
vacancyflux_cahnhilliard_thermalFluc vacancyflux_cahnhilliard_thermalFluc
real(pReal), parameter, private :: & real(pReal), parameter, private :: &