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

Merge remote-tracking branch 'origin/development' into development

This commit is contained in:
Philip Eisenlohr 2017-07-31 09:51:10 -04:00
parent 837b14592a a4ceebc00b
commit c8222d6f82
24 changed files with 382 additions and 273 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
.gitlab-ci.yml
View File

@ -1,3 +1,4 @@
---
stages:
- prepareAll
- preprocessing
@ -23,26 +24,30 @@ stages:
###################################################################################################
before_script:
- if [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue)x == 'x' ]; then echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue; fi
- while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) -ne 1 ];do sleep 5m; done
- source $DAMASKROOT/DAMASK_env.sh
- if [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue)x == 'x' ];
then echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue;
fi
- while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ];
do sleep 5m;
done
- source $DAMASKROOT/env/DAMASK.sh
- cd $DAMASKROOT/PRIVATE/testing
###################################################################################################
variables:
#================================================================================================
# ===============================================================================================
# GitLab Settings
#================================================================================================
# ===============================================================================================
GIT_SUBMODULE_STRATEGY: none
#================================================================================================
# ===============================================================================================
# Shortcut names
#================================================================================================
DAMASKROOT: "$TESTROOT/GitLabCI_Pipeline_$CI_PIPELINE_ID/DAMASK"
# ===============================================================================================
DAMASKROOT: "$TESTROOT/GitLabCI_Pipeline_$CI_PIPELINE_ID/DAMASK"
#================================================================================================
# ===============================================================================================
# Names of module files to load
#================================================================================================
# ===============================================================================================
# ++++++++++++ Compiler ++++++++++++++++++++++++++++++++++++++++++++++
IntelCompiler16_0: "Compiler/Intel/16.0 Libraries/IMKL/2016"
IntelCompiler17_0: "Compiler/Intel/17.0 Libraries/IMKL/2017"
@ -85,14 +90,17 @@ checkout:
stage: prepareAll
before_script:
- echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue
- while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) -ne 1 ];do sleep 5m; done
- while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ];
do sleep 5m;
done
script:
- mkdir -p $DAMASKROOT
- cd $DAMASKROOT
- if [ -d DAMASK ]; then rm -rf DAMASK; fi # there might be some leftovers from a failed clone
- git clone -q git@magit1.mpie.de:damask/DAMASK.git .
- git checkout $CI_COMMIT_SHA
- git submodule update --init
- source DAMASK_env.sh
- source env/DAMASK.sh
- make processing
except:
- master

2
VERSION
View File

@ -1 +1 @@
v2.0.1-728-g66efd90
v2.0.1-805-gdc3eda3

21
env/DAMASK.csh vendored
View File

@ -3,16 +3,21 @@
set CALLED=($_)
set DIRNAME=`dirname $CALLED[2]`
# transition compatibility (renamed $DAMASK_ROOT/DAMASK_env.csh to $DAMASK_ROOT/env/DAMASK.csh)
set FILENAME=`basename $CALLED[2]`
if ($FILENAME == "DAMASK.csh") then
set DIRNAME=$DIRNAME"/../"
endif
set DAMASK_ROOT=`python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" $DIRNAME`
source $DAMASK_ROOT/CONFIG
# if DAMASK_BIN is present and not in $PATH, add it
# if DAMASK_BIN is present
if ( $?DAMASK_BIN) then
set MATCH=`echo :${PATH}: | grep ${DAMASK_BIN}:`
if ( "x$MATCH" == "x" ) then
set PATH=${DAMASK_BIN}:${PATH}
endif
set path = ($DAMASK_BIN $path)
endif
set SOLVER=`which DAMASK_spectral`
@ -62,4 +67,8 @@ if ( $?prompt ) then
endif
setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS
setenv PYTHONPATH $DAMASK_ROOT/lib:$PYTHONPATH
if ( ! $?PYTHONPATH ) then
setenv PYTHONPATH $DAMASK_ROOT/lib
else
setenv PYTHONPATH $DAMASK_ROOT/lib:$PYTHONPATH
endif

29
env/DAMASK.sh vendored
View File

@ -1,5 +1,5 @@
# sets up an environment for DAMASK on bash
# usage: source DAMASK_env.sh
# usage: source DAMASK.sh
if [ "$OSTYPE" == "linux-gnu" ] || [ "$OSTYPE" == 'linux' ]; then
@ -10,8 +10,13 @@ else
DAMASK_ROOT=${STAT##* }
fi
# transition compatibility (renamed $DAMASK_ROOT/DAMASK_env.sh to $DAMASK_ROOT/env/DAMASK.sh)
if [ ${BASH_SOURCE##*/} == "DAMASK.sh" ]; then
DAMASK_ROOT=$DAMASK_ROOT'/..'
fi
# shorthand command to change to DAMASK_ROOT directory
eval "function damask() { cd $DAMASK_ROOT; }"
eval "function DAMASK_root() { cd $DAMASK_ROOT; }"
# defining set() allows to source the same file for tcsh and bash, with and without space around =
set() {
@ -20,22 +25,16 @@ set() {
source $DAMASK_ROOT/CONFIG
unset -f set
# add DAMASK_BIN if present but not yet in $PATH
if [[ "x$DAMASK_BIN" != "x" && ! $(echo ":$PATH:" | grep $DAMASK_BIN:) ]]; then
export PATH=$DAMASK_BIN:$PATH
fi
# add DAMASK_BIN if present
[ "x$DAMASK_BIN" != "x" ] && PATH=$DAMASK_BIN:$PATH
SOLVER=$(which DAMASK_spectral || true 2>/dev/null)
if [ "x$SOLVER" == "x" ]; then
SOLVER='Not found!'
fi
[ "x$SOLVER" == "x" ] && SOLVER='Not found!'
PROCESSING=$(which postResults || true 2>/dev/null)
if [ "x$PROCESSING" == "x" ]; then
PROCESSING='Not found!'
fi
if [ "x$DAMASK_NUM_THREADS" == "x" ]; then
DAMASK_NUM_THREADS=1
fi
[ "x$PROCESSING" == "x" ] && PROCESSING='Not found!'
[ "x$DAMASK_NUM_THREADS" == "x" ] && DAMASK_NUM_THREADS=1
# according to http://software.intel.com/en-us/forums/topic/501500
# this seems to make sense for the stack size

22
env/DAMASK.zsh vendored
View File

@ -1,10 +1,15 @@
# sets up an environment for DAMASK on zsh
# usage: source DAMASK_env.zsh
# usage: source DAMASK.zsh
DAMASK_ROOT=${0:a:h}
# transition compatibility (renamed $DAMASK_ROOT/DAMASK_env.zsh to $DAMASK_ROOT/env/DAMASK.zsh)
if [ ${0:t:r} = 'DAMASK' ]; then
DAMASK_ROOT=${0:a:h}'/..'
else
DAMASK_ROOT=${0:a:h}
fi
# shorthand command to change to DAMASK_ROOT directory
eval "function damask() { cd $DAMASK_ROOT; }"
eval "function DAMASK_root() { cd $DAMASK_ROOT; }"
# defining set() allows to source the same file for tcsh and zsh, with and without space around =
set() {
@ -13,17 +18,12 @@ set() {
source $DAMASK_ROOT/CONFIG
unset -f set
# add DAMASK_BIN if present but not yet in $PATH
MATCH=`echo ":$PATH:" | grep $DAMASK_BIN:`
if [[ ( "x$DAMASK_BIN" != "x" ) && ( "x$MATCH" = "x" ) ]]; then
export PATH=$DAMASK_BIN:$PATH
fi
# add DAMASK_BIN if present
[ "x$DAMASK_BIN" != "x" ] && PATH=$DAMASK_BIN:$PATH
SOLVER=`which DAMASK_spectral || True 2>/dev/null`
PROCESSING=`which postResults || True 2>/dev/null`
if [ "x$DAMASK_NUM_THREADS" = "x" ]; then
DAMASK_NUM_THREADS=1
fi
[ "x$DAMASK_NUM_THREADS" = "x" ] && DAMASK_NUM_THREADS=1
# according to http://software.intel.com/en-us/forums/topic/501500
# this seems to make sense for the stack size

2
examples/ConfigFiles/Kinematics_Thermal_Expansion.config
View File

@ -1,2 +1,2 @@
(kinematics) thermal_expansion
thermal_expansion11 0.00231
thermal_expansion11 0.0000231

2
examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config
View File

@ -43,7 +43,7 @@ q_slip 1.0 # q-exponent in glide velocity
CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path
D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s]
Qsd 4.5e-19 # Activation energy for climb [J]
Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b]
Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b^3]
Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b]
atol_rho 1.0
interactionSlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Interaction coefficients (Kubin et al. 2008)

6
installation/mods_MarcMentat/apply_DAMASK_modifications.sh
View File

@ -102,8 +102,7 @@ for filename in 'edit_window' \
done
# Mentat scripts
echo ''
echo 'adapting Mentat menus...'
echo -e '\nadapting Mentat menus...'
theDIR=$INSTALLDIR/mentat$VERSION/menus
for filename in 'job_run.ms'; do
cp $SCRIPTLOCATION/$VERSION/Mentat_menus/$filename $theDIR
@ -163,5 +162,4 @@ echo ''
echo 'precompiling $VERSION HYPELA2 user subroutine...'
echo 'not yet implemented..!'
echo ''
echo 'done.'
echo -e '\ndone.'

2
installation/symlink_Processing.py
View File

@ -30,7 +30,7 @@ for subDir in processing_subDirs:
for theFile in os.listdir(theDir):
theName,theExt = os.path.splitext(theFile)
if theExt in processing_extensions: # only consider files with proper extensions
if theExt in processing_extensions: # only consider files with proper extensions
src = os.path.abspath(os.path.join(theDir,theFile))
sym_link = os.path.abspath(os.path.join(binDir,theName))

137
processing/post/addGaussian.py Executable file
View File

@ -0,0 +1,137 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys
import numpy as np
from optparse import OptionParser
from scipy import ndimage
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
Add column(s) containing Gaussian filtered values of requested column(s).
Operates on periodic and non-periodic ordered three-dimensional data sets.
For Details see scipy.ndimage documentation.
""", version = scriptID)
parser.add_option('-p','--pos','--periodiccellcenter',
dest = 'pos',
type = 'string', metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('-s','--scalar',
dest = 'scalar',
action = 'extend', metavar = '<string LIST>',
help = 'label(s) of scalar field values')
parser.add_option('-o','--order',
dest = 'order',
type = int,
metavar = 'int',
help = 'order of the filter')
parser.add_option('--sigma',
dest = 'sigma',
type = float,
metavar = 'float',
help = 'standard deviation')
parser.add_option('--periodic',
dest = 'periodic',
action = 'store_true',
help = 'assume periodic grain structure'
)
parser.set_defaults(pos = 'pos',
order = 0,
sigma = 1,
periodic = False
)
(options,filenames) = parser.parse_args()
if options.scalar is None:
parser.error('no data column specified.')
# --- loop over input files ------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try: table = damask.ASCIItable(name = name,buffered = False)
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
items = {
'scalar': {'dim': 1, 'shape': [1], 'labels':options.scalar, 'active':[], 'column': []},
}
errors = []
remarks = []
column = {}
if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos))
else: colCoord = table.label_index(options.pos)
for type, data in items.iteritems():
for what in (data['labels'] if data['labels'] is not None else []):
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
else:
items[type]['active'].append(what)
items[type]['column'].append(table.label_index(what))
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
for type, data in items.iteritems():
for label in data['active']:
table.labels_append(['Gauss{}({})'.format(options.sigma,label)]) # extend ASCII header with new labels
table.head_write()
# --------------- figure out size and grid ---------------------------------------------------------
table.data_readArray()
coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones
# ------------------------------------------ process value field -----------------------------------
stack = [table.data]
for type, data in items.iteritems():
for i,label in enumerate(data['active']):
stack.append(ndimage.filters.gaussian_filter(table.data[:,data['column'][i]],
options.sigma,options.order,
mode = 'wrap' if options.periodic else 'nearest'
).reshape([table.data.shape[0],1])
)
# ------------------------------------------ output result -----------------------------------------
if len(stack) > 1: table.data = np.hstack(tuple(stack))
table.data_writeArray('%.12g')
# ------------------------------------------ output finalization -----------------------------------
table.close() # close input ASCII table (works for stdin)

2
processing/post/addGrainID.py
View File

@ -148,7 +148,7 @@ for name in filenames:
bg.set_message('reading positions...')
table.data_readArray(options.pos) # read position vectors
table.data_readArray(options.pos) # read position vectors
grainID = -np.ones(len(table.data),dtype=int)
start = tick = time.clock()

3
processing/pre/geom_fromTable.py
View File

@ -23,7 +23,7 @@ Generate geometry description and material configuration from position, phase, a
parser.add_option('--coordinates',
dest = 'pos',
type = 'string', metavar = 'string',
help = 'coordinates label')
help = 'coordinates label (%default)')
parser.add_option('--phase',
dest = 'phase',
type = 'string', metavar = 'string',
@ -90,6 +90,7 @@ parser.set_defaults(symmetry = [damask.Symmetry.lattices[-1]],
homogenization = 1,
crystallite = 1,
verbose = False,
pos = 'pos',
)
(options,filenames) = parser.parse_args()

0
processing/pre/mentat_pbcOnBoxMesh.py Normal file → Executable file
View File

15
src/CMakeLists.txt
View File

@ -1,16 +1,5 @@
# special flags for some files
if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU")
SET_SOURCE_FILES_PROPERTIES( "prec.f90" PROPERTIES
COMPILE_FLAGS "-fno-range-check -fall-intrinsics -fno-fast-math")
# fno-range-check: Disable range checking on results of simplification of constant expressions during compilation
# --> allows the definition of DAMASK_NaN
#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored
# and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external
# --> allows the use of 'isnan'
#-fno-fast-math:
# --> otherwise, when setting -ffast-math, isnan always evaluates to false (I would call it a bug)
SET_SOURCE_FILES_PROPERTIES( "lattice.f90" PROPERTIES
COMPILE_FLAGS "-ffree-line-length-240")
# long lines for interaction matrix
@ -185,8 +174,10 @@ if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral")
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>)
if(NOT "${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY")
add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES})
add_dependencies(DAMASK_spectral SPECTRAL_SOLVER)
else()
add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90")
endif()
add_dependencies(DAMASK_spectral SPECTRAL_SOLVER)
elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM")
add_library(FEM_UTILITIES OBJECT "FEM_utilities.f90")
add_dependencies(FEM_UTILITIES DAMASK_CPFE)

138
src/DAMASK_spectral.f90 Normal file → Executable file
View File

@ -364,12 +364,12 @@ program DAMASK_spectral
case (DAMASK_spectral_SolverBasicPETSc_label)
call basicPETSc_init
case (DAMASK_spectral_SolverAL_label)
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) &
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
call AL_init
case (DAMASK_spectral_SolverPolarisation_label)
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) &
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
call Polarisation_init
@ -450,8 +450,7 @@ program DAMASK_spectral
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
if (worldrank == 0) &
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
endif
!--------------------------------------------------------------------------------------------------
! loopping over loadcases
@ -498,22 +497,20 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! 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
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
!--------------------------------------------------------------------------------------------------
! forward fields
@ -595,7 +592,7 @@ program DAMASK_spectral
cutBack = .False.
if(solres(1)%termIll .or. .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'
write(6,'(/,a)') ' cut back detected'
cutBack = .True.
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt
@ -603,11 +600,15 @@ program DAMASK_spectral
timeinc = timeinc/2.0_pReal
elseif (solres(1)%termIll) then ! material point model cannot find a solution, exit in any casy
call IO_warning(850_pInt)
call MPI_file_close(resUnit,ierr)
close(statUnit)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
elseif (continueCalculation == 1_pInt) then
guess = .true. ! accept non converged BVP solution
else ! default behavior, exit if spectral solver does not converge
call IO_warning(850_pInt)
call MPI_file_close(resUnit,ierr)
close(statUnit)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
endif
else
@ -624,13 +625,11 @@ program DAMASK_spectral
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) &
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &
' increment ', totalIncsCounter, ' converged'
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &
' increment ', totalIncsCounter, ' converged'
else
if (worldrank == 0) &
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
@ -665,45 +664,15 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! 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
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!'
call MPI_file_close(resUnit,ierr)
close(statUnit)
do field = 1, nActiveFields
select case(loadCases(1)%ID(field))
case(FIELD_MECH_ID)
select case (spectral_solver)
case (DAMASK_spectral_SolverBasicPETSc_label)
call BasicPETSC_destroy()
case (DAMASK_spectral_SolverAL_label)
call AL_destroy()
case (DAMASK_spectral_SolverPolarisation_label)
call Polarisation_destroy()
end select
case(FIELD_THERMAL_ID)
call spectral_thermal_destroy()
case(FIELD_DAMAGE_ID)
call spectral_damage_destroy()
end select
enddo
call utilities_destroy()
call PETScFinalize(ierr); CHKERRQ(ierr)
#ifdef _OPENMP
call MPI_finalize(i)
if (i /= 0_pInt) then
call IO_error(error_ID=894, el=i, ext_msg="Finalize()")
endif
#endif
if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;)
@ -721,10 +690,46 @@ end program DAMASK_spectral
subroutine quit(stop_id)
use prec, only: &
pInt
use spectral_mech_Basic, only: &
BasicPETSC_destroy
use spectral_mech_AL, only: &
AL_destroy
use spectral_mech_Polarisation, only: &
Polarisation_destroy
use spectral_damage, only: &
spectral_damage_destroy
use spectral_thermal, only: &
spectral_thermal_destroy
use spectral_utilities, only: &
utilities_destroy
implicit none
#include <petsc/finclude/petscsys.h>
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
integer(pInt) :: error = 0_pInt
PetscErrorCode :: ierr = 0
logical :: ErrorInQuit
external :: &
PETScFinalize, &
MPI_finalize
call BasicPETSC_destroy()
call AL_destroy()
call Polarisation_destroy()
call spectral_damage_destroy()
call spectral_thermal_destroy()
call utilities_destroy()
call PETScFinalize(ierr)
if(ierr /= 0) write(6,'(a)') ' Error in PETScFinalize'
#ifdef _OPENMP
call MPI_finalize(error)
if(error /= 0) write(6,'(a)') ' Error in MPI_finalize'
#endif
ErrorInQuit = (ierr /= 0 .or. error /= 0_pInt)
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
@ -735,12 +740,13 @@ subroutine quit(stop_id)
dateAndTime(6),':',&
dateAndTime(7)
if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id < 0_pInt) then ! terminally ill, restart might help
if (stop_id == 0_pInt .and. .not. ErrorInQuit) stop 0 ! normal termination
if (stop_id < 0_pInt .and. .not. ErrorInQuit) then ! terminally ill, restart might help
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
if (stop_id == 3_pInt .and. .not. ErrorInQuit) stop 3 ! not all incs converged
stop 1 ! error (message from IO_error)
end subroutine quit

90
src/crystallite.f90
View File

@ -1232,8 +1232,8 @@ end subroutine crystallite_stressAndItsTangent
!> @brief integrate stress, state with 4th order explicit Runge Kutta method
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRK4()
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use numerics, only: &
numerics_integrationMode
use debug, only: &
@ -1331,9 +1331,9 @@ subroutine crystallite_integrateStateRK4()
if (crystallite_todo(g,i,e)) then
c = phasememberAt(g,i,e)
p = phaseAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -1475,9 +1475,9 @@ subroutine crystallite_integrateStateRK4()
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -1528,8 +1528,8 @@ end subroutine crystallite_integrateStateRK4
!> adaptive step size (use 5th order solution to advance = "local extrapolation")
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRKCK45()
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
debug_level, &
debug_crystallite, &
@ -1647,9 +1647,9 @@ subroutine crystallite_integrateStateRKCK45()
if (crystallite_todo(g,i,e)) then
cc = phasememberAt(g,i,e)
p = phaseAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -1801,9 +1801,9 @@ subroutine crystallite_integrateStateRKCK45()
p = phaseAt(g,i,e)
cc = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -2031,8 +2031,8 @@ end subroutine crystallite_integrateStateRKCK45
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateAdaptiveEuler()
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
debug_level, &
debug_crystallite, &
@ -2133,9 +2133,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -2254,9 +2254,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -2391,8 +2391,8 @@ end subroutine crystallite_integrateStateAdaptiveEuler
!> @brief integrate stress, and state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateEuler()
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
debug_level, &
debug_crystallite, &
@ -2471,9 +2471,9 @@ eIter = FEsolving_execElem(1:2)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
c = phasememberAt(g,i,e)
p = phaseAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local...
@ -2614,8 +2614,8 @@ end subroutine crystallite_integrateStateEuler
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateFPI()
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
debug_e, &
debug_i, &
@ -2737,9 +2737,9 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
@ -2851,9 +2851,9 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
crystallite_todo(g,i,e) = .false. ! ... skip me next time
@ -3062,8 +3062,9 @@ end subroutine crystallite_integrateStateFPI
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!--------------------------------------------------------------------------------------------------
logical function crystallite_stateJump(ipc,ip,el)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: &
prec_isNaN, &
dNeq0
use debug, only: &
debug_level, &
@ -3098,7 +3099,7 @@ logical function crystallite_stateJump(ipc,ip,el)
p = phaseAt(ipc,ip,el)
call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe(1:3,1:3,ipc,ip,el), ipc,ip,el)
mySizePlasticDeltaState = plasticState(p)%sizeDeltaState
if( any(prec_isNaN(plasticState(p)%deltaState(:,c)))) then ! NaN occured in deltaState
if( any(IEEE_is_NaN(plasticState(p)%deltaState(:,c)))) then ! NaN occured in deltaState
crystallite_stateJump = .false.
return
endif
@ -3106,7 +3107,7 @@ logical function crystallite_stateJump(ipc,ip,el)
plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
do mySource = 1_pInt, phase_Nsources(p)
mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState
if( any(prec_isNaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState
if( any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState
crystallite_stateJump = .false.
return
endif
@ -3169,15 +3170,18 @@ logical function crystallite_integrateStress(&
el,& ! element number
timeFraction &
)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: pLongInt, &
tol_math_check, &
prec_isNaN, &
dEq0
use numerics, only: nStress, &
aTol_crystalliteStress, &
rTol_crystalliteStress, &
iJacoLpresiduum, &
numerics_integrationMode
numerics_integrationMode, &
subStepSizeLp, &
subStepSizeLi
use debug, only: debug_level, &
debug_crystallite, &
debug_levelBasic, &
@ -3263,9 +3267,7 @@ logical function crystallite_integrateStress(&
dLp_dT3333, &
dLi_dT3333
real(pReal) detInvFi, & ! determinant of InvFi
steplengthLp0, &
steplengthLp, &
steplengthLi0, &
steplengthLi, &
dt, & ! time increment
aTolLp, &
@ -3351,8 +3353,7 @@ logical function crystallite_integrateStress(&
NiterationStressLi = 0_pInt
jacoCounterLi = 0_pInt
steplengthLi0 = 1.0_pReal
steplengthLi = steplengthLi0
steplengthLi = 1.0_pReal
residuumLi_old = 0.0_pReal
LiLoop: do
@ -3372,8 +3373,7 @@ logical function crystallite_integrateStress(&
NiterationStressLp = 0_pInt
jacoCounterLp = 0_pInt
steplengthLp0 = 1.0_pReal
steplengthLp = steplengthLp0
steplengthLp = 1.0_pReal
residuumLp_old = 0.0_pReal
Lpguess_old = Lpguess
@ -3430,7 +3430,7 @@ logical function crystallite_integrateStress(&
aTol_crystalliteStress) ! minimum lower cutoff
residuumLp = Lpguess - Lp_constitutive
if (any(prec_isNaN(residuumLp))) then ! NaN in residuum...
if (any(IEEE_is_NaN(residuumLp))) then ! NaN in residuum...
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip ipc ', &
@ -3445,9 +3445,9 @@ logical function crystallite_integrateStress(&
.or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLp_old = residuumLp ! ...remember old values and...
Lpguess_old = Lpguess
steplengthLp = steplengthLp0 ! ...proceed with normal step length (calculate new search direction)
steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
steplengthLp = 0.5_pReal * steplengthLp ! ...try with smaller step length in same direction
steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplengthLp * deltaLp
cycle LpLoop
endif
@ -3520,7 +3520,7 @@ logical function crystallite_integrateStress(&
aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff
residuumLi = Liguess - Li_constitutive
if (any(prec_isNaN(residuumLi))) then ! NaN in residuum...
if (any(IEEE_is_NaN(residuumLi))) then ! NaN in residuum...
return ! ...me = .false. to inform integrator about problem
elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance
exit LiLoop ! ...leave iteration loop
@ -3528,9 +3528,9 @@ logical function crystallite_integrateStress(&
.or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLi_old = residuumLi ! ...remember old values and...
Liguess_old = Liguess
steplengthLi = steplengthLi0 ! ...proceed with normal step length (calculate new search direction)
steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
steplengthLi = 0.5_pReal * steplengthLi ! ...try with smaller step length in same direction
steplengthLi = subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction
Liguess = Liguess_old + steplengthLi * deltaLi
cycle LiLoop
endif

5
src/homogenization.f90
View File

@ -923,12 +923,11 @@ subroutine materialpoint_postResults
material_phase, &
homogenization_Ngrains, &
microstructure_crystallite
use constitutive, only: &
#ifdef FEM
use constitutive, only: &
constitutive_plasticity_maxSizePostResults, &
constitutive_source_maxSizePostResults, &
constitutive_source_maxSizePostResults
#endif
constitutive_postResults
use crystallite, only: &
#ifdef FEM
crystallite_maxSizePostResults, &

6
src/lattice.f90
View File

@ -2170,8 +2170,8 @@ pure function lattice_symmetrize33(struct,T33)
!> @brief figures whether unit quat falls into stereographic standard triangle
!--------------------------------------------------------------------------------------------------
logical pure function lattice_qInSST(Q, struct)
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use math, only: &
math_qToRodrig
@ -2181,7 +2181,7 @@ logical pure function lattice_qInSST(Q, struct)
real(pReal), dimension(3) :: Rodrig ! Rodrigues vector of Q
Rodrig = math_qToRodrig(Q)
if (any(prec_isNaN(Rodrig))) then
if (any(IEEE_is_NaN(Rodrig))) then
lattice_qInSST = .false.
else
select case (struct)

36
src/math.f90
View File

@ -1666,15 +1666,16 @@ end function math_qToEulerAxisAngle
!> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz)
!--------------------------------------------------------------------------------------------------
pure function math_qToRodrig(Q)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: &
DAMASK_NaN, &
tol_math_check
implicit none
real(pReal), dimension(4), intent(in) :: Q
real(pReal), dimension(3) :: math_qToRodrig
math_qToRodrig = merge(Q(2:4)/Q(1),DAMASK_NaN,abs(Q(1)) > tol_math_check) ! NaN for 180 deg since Rodrig is unbound
math_qToRodrig = merge(Q(2:4)/Q(1),IEEE_value(1.0_pReal,IEEE_quiet_NaN),abs(Q(1)) > tol_math_check)! NaN for 180 deg since Rodrig is unbound
end function math_qToRodrig
@ -2092,11 +2093,11 @@ end function math_rotationalPart33
!--------------------------------------------------------------------------------------------------
!> @brief Eigenvalues of symmetric matrix m
! will return NaN on error
! will return NaN on error
!--------------------------------------------------------------------------------------------------
function math_eigenvaluesSym(m)
use prec, only: &
DAMASK_NaN
use, intrinsic :: &
IEEE_arithmetic
implicit none
real(pReal), dimension(:,:), intent(in) :: m
@ -2109,7 +2110,7 @@ function math_eigenvaluesSym(m)
vectors = m ! copy matrix to input (doubles as output) array
call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info)
if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN
if (info /= 0_pInt) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
end function math_eigenvaluesSym
@ -2701,29 +2702,13 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
end function math_rotate_forward3333
!--------------------------------------------------------------------------------------------------
!> @brief calculate average of tensor field
!--------------------------------------------------------------------------------------------------
function math_tensorAvg(field)
implicit none
real(pReal), dimension(3,3) :: math_tensorAvg
real(pReal), intent(in), dimension(:,:,:,:,:) :: field
real(pReal) :: wgt
wgt = 1.0_pReal/real(size(field,3)*size(field,4)*size(field,5), pReal)
math_tensorAvg = sum(sum(sum(field,dim=5),dim=4),dim=3)*wgt
end function math_tensorAvg
!--------------------------------------------------------------------------------------------------
!> @brief limits a scalar value to a certain range (either one or two sided)
! Will return NaN if left > right
!--------------------------------------------------------------------------------------------------
real(pReal) pure function math_limit(a, left, right)
use prec, only: &
DAMASK_NaN
use, intrinsic :: &
IEEE_arithmetic
implicit none
real(pReal), intent(in) :: a
@ -2735,7 +2720,8 @@ real(pReal) pure function math_limit(a, left, right)
merge(right, huge(a), present(right)) &
)
if (present(left) .and. present(right)) math_limit = merge (DAMASK_NaN,math_limit, left>right)
if (present(left) .and. present(right)) &
math_limit = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_limit, left>right)
end function math_limit

10
src/numerics.f90
View File

@ -42,6 +42,8 @@ module numerics
subStepMinHomog = 1.0e-3_pReal, & !< minimum (relative) size of sub-step allowed during cutback in homogenization
subStepSizeCryst = 0.25_pReal, & !< size of first substep when cutback in crystallite
subStepSizeHomog = 0.25_pReal, & !< size of first substep when cutback in homogenization
subStepSizeLp = 0.5_pReal, & !< size of first substep when cutback in Lp calculation
subStepSizeLi = 0.5_pReal, & !< size of first substep when cutback in Li calculation
stepIncreaseCryst = 1.5_pReal, & !< increase of next substep size when previous substep converged in crystallite
stepIncreaseHomog = 1.5_pReal, & !< increase of next substep size when previous substep converged in homogenization
rTol_crystalliteState = 1.0e-6_pReal, & !< relative tolerance in crystallite state loop
@ -295,6 +297,10 @@ subroutine numerics_init
subStepSizeCryst = IO_floatValue(line,chunkPos,2_pInt)
case ('stepincreasecryst')
stepIncreaseCryst = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizelp')
subStepSizeLp = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizeli')
subStepSizeLi = IO_floatValue(line,chunkPos,2_pInt)
case ('substepminhomog')
subStepMinHomog = IO_floatValue(line,chunkPos,2_pInt)
case ('substepsizehomog')
@ -515,6 +521,8 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' subStepMinCryst: ',subStepMinCryst
write(6,'(a24,1x,es8.1)') ' subStepSizeCryst: ',subStepSizeCryst
write(6,'(a24,1x,es8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst
write(6,'(a24,1x,es8.1)') ' subStepSizeLp: ',subStepSizeLp
write(6,'(a24,1x,es8.1)') ' subStepSizeLi: ',subStepSizeLi
write(6,'(a24,1x,i8)') ' nState: ',nState
write(6,'(a24,1x,i8)') ' nStress: ',nStress
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
@ -643,6 +651,8 @@ subroutine numerics_init
if (subStepMinCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinCryst')
if (subStepSizeCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeCryst')
if (stepIncreaseCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseCryst')
if (subStepSizeLp <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLp')
if (subStepSizeLi <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLi')
if (subStepMinHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinHomog')
if (subStepSizeHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeHomog')
if (stepIncreaseHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseHomog')

17
src/plastic_isotropic.f90
View File

@ -9,8 +9,7 @@
module plastic_isotropic
use prec, only: &
pReal,&
pInt, &
DAMASK_NaN
pInt
implicit none
private
@ -36,14 +35,14 @@ module plastic_isotropic
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
real(pReal) :: &
fTaylor = DAMASK_NaN, &
tau0 = DAMASK_NaN, &
gdot0 = DAMASK_NaN, &
n = DAMASK_NaN, &
h0 = DAMASK_NaN, &
fTaylor, &
tau0, &
gdot0, &
n, &
h0, &
h0_slopeLnRate = 0.0_pReal, &
tausat = DAMASK_NaN, &
a = DAMASK_NaN, &
tausat, &
a, &
aTolFlowstress = 1.0_pReal, &
aTolShear = 1.0e-6_pReal, &
tausat_SinhFitA= 0.0_pReal, &

10
src/plastic_nonlocal.f90
View File

@ -2377,9 +2377,9 @@ end subroutine plastic_nonlocal_deltaState
!---------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
timestep,subfrac, ip,el)
use prec, only: DAMASK_NaN, &
dNeq0, &
use, intrinsic :: &
IEEE_arithmetic
use prec, only: dNeq0, &
dNeq, &
dEq0
use numerics, only: numerics_integrationMode, &
@ -2701,7 +2701,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
write(6,'(a)') '<< CONST >> enforcing cutback !!!'
endif
#endif
plasticState(p)%dotState = DAMASK_NaN ! -> return NaN and, hence, enforce cutback
plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback
return
endif
@ -2984,7 +2984,7 @@ if ( any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(in
write(6,'(a)') '<< CONST >> enforcing cutback !!!'
endif
#endif
plasticState(p)%dotState = DAMASK_NaN
plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
return
else
forall (s = 1:ns, t = 1_pInt:4_pInt)

36
src/prec.f90
View File

@ -5,22 +5,12 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @brief setting precision for real and int type
!> @details setting precision for real and int type and for DAMASK_NaN. Definition is made
!! depending on makro "INT" defined during compilation
!! for details on NaN see https://software.intel.com/en-us/forums/topic/294680
!--------------------------------------------------------------------------------------------------
module prec
#if !(defined(__GFORTRAN__) && __GNUC__ < 5)
use, intrinsic :: & ! unfortunately not avialable in gfortran <= 5
IEEE_arithmetic
#endif
implicit none
private
#if (FLOAT==8)
integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300)
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000',pReal) !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html)
#else
NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION
#endif
@ -106,7 +96,6 @@ module prec
public :: &
prec_init, &
prec_isNaN, &
dEq, &
dEq0, &
cEq, &
@ -118,7 +107,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief reporting precision and checking if DAMASK_NaN is set correctly
!> @brief reporting precision
!--------------------------------------------------------------------------------------------------
subroutine prec_init
use, intrinsic :: &
@ -133,36 +122,13 @@ subroutine prec_init
write(6,'(a,i3)') ' Bytes for pReal: ',pReal
write(6,'(a,i3)') ' Bytes for pInt: ',pInt
write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt
write(6,'(a,e10.3)') ' NaN: ', DAMASK_NaN
write(6,'(a,l3)') ' NaN != NaN: ',DAMASK_NaN /= DAMASK_NaN
write(6,'(a,l3,/)') ' NaN check passed ',prec_isNAN(DAMASK_NaN)
if ((.not. prec_isNaN(DAMASK_NaN)) .or. (DAMASK_NaN == DAMASK_NaN)) call quit(9000)
realloc_lhs_test = [1_pInt,2_pInt]
if (realloc_lhs_test(2)/=2_pInt) call quit(9000)
end subroutine prec_init
!--------------------------------------------------------------------------------------------------
!> @brief figures out if a floating point number is NaN
! basically just a small wrapper, because gfortran < 5.0 does not have the IEEE module
!--------------------------------------------------------------------------------------------------
logical elemental pure function prec_isNaN(a)
implicit none
real(pReal), intent(in) :: a
#if (defined(__GFORTRAN__) && __GNUC__ < 5)
intrinsic :: isNaN
prec_isNaN = isNaN(a)
#else
prec_isNaN = IEEE_is_NaN(a)
#endif
end function prec_isNaN
!--------------------------------------------------------------------------------------------------
!> @brief equality comparison for float with double precision
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq

6
src/spectral_utilities.f90
View File

@ -741,8 +741,8 @@ end function utilities_curlRMS
!> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC
!--------------------------------------------------------------------------------------------------
function utilities_maskedCompliance(rot_BC,mask_stress,C)
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use IO, only: &
IO_error
use math, only: &
@ -794,7 +794,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
if (any(prec_isNaN(s_reduced))) errmatinv = .true.
if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
temp99_Real = 0.0_pReal ! fill up compliance with zeros
k = 0_pInt