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 MiscImprovements

This commit is contained in:
Martin Diehl 2020-02-11 17:11:43 +01:00
parent cd76172d74 3dacd3d931
commit c2c84d698f
25 changed files with 496 additions and 4098 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
.gitlab-ci.yml
View File

@ -75,7 +75,6 @@ variables:
# ------------ Defaults ---------------------------------------------- # ------------ Defaults ----------------------------------------------
MSC: "$MSC2019" MSC: "$MSC2019"
IntelMarc: "$IntelCompiler17_8" IntelMarc: "$IntelCompiler17_8"
IntelAbaqus: "$IntelCompiler16_4"
HDF5Marc: "HDF5/1.10.5/Intel-17.8" HDF5Marc: "HDF5/1.10.5/Intel-17.8"
# ++++++++++++ Documentation ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++ Documentation ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Doxygen1_8_15: "Documentation/Doxygen/1.8.15" Doxygen1_8_15: "Documentation/Doxygen/1.8.15"
@ -476,15 +475,6 @@ createTar:
- release - release
################################################################################################### ###################################################################################################
AbaqusStd:
stage: createDocumentation
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
- $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT abaqus
except:
- master
- release
Marc: Marc:
stage: createDocumentation stage: createDocumentation
script: script:

6
DAMASK_prerequisites.sh
View File

@ -122,9 +122,3 @@ done
firstLevel "CMake" firstLevel "CMake"
getDetails cmake --version getDetails cmake --version
firstLevel "Abaqus"
cd installation/mods_Abaqus # to have the right environment file
for executable in abaqus abq2017 abq2018 abq2019; do
getDetails $executable 'information=all'
done
cd ../..

2
PRIVATE

@ -1 +1 @@
Subproject commit 64432754ce3c590c882cf4987695539cee524ee8 Subproject commit 9dc7065beedf2a097aa60a656bbfa52e55d7147c

2
VERSION
View File

@ -1 +1 @@
v2.0.3-1624-g47109b90 v2.0.3-1639-g81ae6686

2
env/CONFIG vendored
View File

@ -3,5 +3,3 @@ set DAMASK_NUM_THREADS = 4
set MSC_ROOT = /opt/msc set MSC_ROOT = /opt/msc
set MARC_VERSION = 2019 set MARC_VERSION = 2019
set ABAQUS_VERSION = 2019

1
env/DAMASK.csh vendored
View File

@ -51,5 +51,4 @@ else
setenv PYTHONPATH $DAMASK_ROOT/python:$PYTHONPATH setenv PYTHONPATH $DAMASK_ROOT/python:$PYTHONPATH
endif endif
setenv MSC_ROOT setenv MSC_ROOT
setenv ABAQUS_VERSION
setenv MARC_VERSION setenv MARC_VERSION

1
examples/ConfigFiles/debug.config
View File

@ -17,7 +17,6 @@ homogenization # homogenization_*.f90 possible values: ba
CPFEM # CPFEM.f90 possible values: basic, extensive, selective CPFEM # CPFEM.f90 possible values: basic, extensive, selective
spectral # DAMASK_spectral.f90 possible values: basic, fft, restart, divergence, rotation, petsc spectral # DAMASK_spectral.f90 possible values: basic, fft, restart, divergence, rotation, petsc
marc # MSC.MARC FEM solver possible values: basic marc # MSC.MARC FEM solver possible values: basic
abaqus # ABAQUS FEM solver possible values: basic
# #
# Parameters for selective # Parameters for selective
element 1 # selected element for debugging (synonymous: "el", "e") element 1 # selected element for debugging (synonymous: "el", "e")

2
examples/ConfigFiles/numerics.config
View File

@ -9,7 +9,7 @@ pert_method 1 # perturbation method (1 = forward, 2 = b
integrator 1 # integration method (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp) integrator 1 # integration method (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp)
integratorStiffness 1 # integration method used for stiffness (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp) integratorStiffness 1 # integration method used for stiffness (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp)
unitlength 1 # physical length of one computational length unit unitlength 1 # physical length of one computational length unit
usepingpong 1 # use the ping pong (collect <-> calc) scheme (always off for Abaqus exp, must be on for Spectral Solver) usepingpong 1 # use the ping pong (collect <-> calc) scheme
## crystallite numerical parameters ## ## crystallite numerical parameters ##
nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")

56
installation/mods_Abaqus/abaqus_v6.env
View File

@ -1,56 +0,0 @@
#
# DAMASK Abaqus Environment File
#
# ------------------------------------
# originally taken from Abaqus ver. 6.11.1
#
#
# Linux (Opteron/EM64T) Settings:
#
# Compile and Link command for user subroutines.
# Compile_cpp and link_exe for Abaqus make utility.
#
import os, re, glob, driverUtils
if False:
# use hdf5 compiler wrapper in $PATH
fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string
link_sl += fortCmd.split()[1:]
fortCmd +=" -DDAMASK_HDF5"
else:
# Use the version in $PATH
fortCmd = "ifort"
# -free to use free-format FORTRAN 90 syntax
# -O <0-3> optimization level
# -fpp use FORTRAN preprocessor on source code
# -fopenmp build with openMP support
# -w90 -w95 suppress messages about use of non-standard Fortran (previous version of abaqus_v6.env only)
# -WB turn a compile-time bounds check into a warning (previous version of abaqus_v6.env only)
# -mP2OPT_hpo_vec_divbyzero=F inofficial compiler switch, proposed by abaqus but highly dubios (previous version of abaqus_v6.env only)
# -ftz flush underflow to zero
# -diag-disable 5268 disable warnings about line length > 132 (only comments there anyway)
# -implicitnone assume no implicit types (e.g. i for integer)
# -standard-semantics sets standard (Fortran 2008) and some other conventions
# -assume nostd_mod_proc_name avoid problems with libraries compiled without that option
# -real-size 64 assume size of real to be 8 bytes, matches our definition of pReal
compile_fortran = (fortCmd + " -c -fPIC -auto -shared-intel " +
"-I%I -free -O3 -fpp -fopenmp " +
"-ftz -diag-disable 5268 " +
"-implicitnone -standard-semantics " +
"-assume nostd_mod_proc_name " +
"-real-size 64 " +
'-DDAMASKVERSION=\\\"n/a\\\"')
# Abaqus/CAE will generate an input file without parts and assemblies.
cae_no_parts_input_file=ON
# Both the Abaqus/Explicit packager and analysis are run in double precision.
double_precision=BOTH
# The user will not be asked whether old job files of the same name should be deleted.
ask_delete=OFF
# usub_lib_dir='your_prefered_location/abqlib'
# Remove the temporary names from the namespace
del fortCmd

61
installation/mods_Abaqus/abaqus_v6_debug.env
View File

@ -1,61 +0,0 @@
#
# DAMASK Abaqus Environment File
#
# ------------------------------------
# originally taken from Abaqus ver. 6.11.1
#
#
# Linux (Opteron/EM64T) Settings:
#
# Compile and Link command for user subroutines.
# Compile_cpp and link_exe for Abaqus make utility.
#
import os, re, glob, driverUtils
if False:
# use hdf5 compiler wrapper in $PATH
fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string
link_sl += fortCmd.split()[1:]
fortCmd +=" -DDAMASK_HDF5"
else:
# Use the version in $PATH
fortCmd = "ifort"
# -free to use free-format FORTRAN 90 syntax
# -O <0-3> optimization level
# -fpp use FORTRAN preprocessor on source code
# -fopenmp build with openMP support
# -w90 -w95 suppress messages about use of non-standard Fortran (previous version of abaqus_v6.env only)
# -WB turn a compile-time bounds check into a warning (previous version of abaqus_v6.env only)
# -mP2OPT_hpo_vec_divbyzero=F inofficial compiler switch, proposed by abaqus but highly dubios (previous version of abaqus_v6.env only)
# -ftz flush underflow to zero
# -diag-disable 5268 disable warnings about line length > 132 (only comments there anyway)
# -implicitnone assume no implicit types (e.g. i for integer)
# -standard-semantics sets standard (Fortran 2008) and some other conventions
# -assume nostd_mod_proc_name avoid problems with libraries compiled without that option
# -real-size 64 assume size of real to be 8 bytes, matches our definition of pReal
# 'check pointers' does not work
compile_fortran = (fortCmd + " -c -fPIC -auto -shared-intel " +
"-I%I -free -O0 -fpp " +
"-ftz -diag-disable 5268 " +
"-implicitnone -standard-semantics " +
"-assume nostd_mod_proc_name " +
"-real-size 64 " +
"-check bounds,format,output_conversion,uninit " +
"-ftrapuv -fpe-all0 " +
"-g -traceback -gen-interfaces -fp-stack-check -fp-model strict " +
'-DDAMASKVERSION=\\\"n/a\\\"')
# Abaqus/CAE will generate an input file without parts and assemblies.
cae_no_parts_input_file=ON
# Both the Abaqus/Explicit packager and analysis are run in double precision.
double_precision=BOTH
# The user will not be asked whether old job files of the same name should be deleted.
ask_delete=OFF
# usub_lib_dir='your_prefered_location/abqlib'
# Remove the temporary names from the namespace
del fortCmd

1
python/damask/environment.py
View File

@ -20,6 +20,5 @@ class Environment():
for item in ['DAMASK_NUM_THREADS', for item in ['DAMASK_NUM_THREADS',
'MSC_ROOT', 'MSC_ROOT',
'MARC_VERSION', 'MARC_VERSION',
'ABAQUS_VERSION',
]: ]:
self.options[item] = os.environ[item] if item in os.environ else None self.options[item] = os.environ[item] if item in os.environ else None

1
python/damask/solver/__init__.py
View File

@ -2,4 +2,3 @@
from .solver import Solver # noqa from .solver import Solver # noqa
from .marc import Marc # noqa from .marc import Marc # noqa
from .abaqus import Abaqus # noqa

36
python/damask/solver/abaqus.py
View File

@ -1,36 +0,0 @@
import subprocess
from .solver import Solver
import damask
class Abaqus(Solver):
"""Wrapper to run DAMASK with Abaqus."""
def __init__(self,version=damask.Environment().options['ABAQUS_VERSION']):
"""
Create a Abaqus solver object.
Parameters
----------
version : integer
Abaqus version
"""
self.solver = 'Abaqus'
try:
self.version = int(version)
except TypeError:
self.version = -1
def return_run_command(self,model):
try:
cmd = 'abq{}'.format(self.version)
subprocess.check_output([cmd,'information=release'])
except OSError: # link to abqXXX not existing
cmd = 'abaqus'
process = subprocess.Popen([cmd,'information=release'],stdout = subprocess.PIPE,stderr = subprocess.PIPE)
detectedVersion = int(process.stdout.readlines()[1].split()[1].decode('utf-8'))
if self.version not in [detectedVersion,-1]:
raise Exception('found Abaqus version {}, but requested {}'.format(detectedVersion,self.version))
return '{} -job {} -user {}/src/DAMASK_abaqus interactive'.format(cmd,model,damask.Environment().rootDir())

35
src/CMakeLists.txt
View File

@ -4,37 +4,36 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
SET_SOURCE_FILES_PROPERTIES("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240") SET_SOURCE_FILES_PROPERTIES("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240")
endif() endif()
file (GLOB damask-sources *.f90 *.c) file(GLOB damask-sources *.f90 *.c)
# probably we should have subfolders for abaqus and MSC.Marc # probably we should have a subfolder for MSC.Marc
list (FILTER damask-sources EXCLUDE REGEX ".*CPFEM.f90") list(FILTER damask-sources EXCLUDE REGEX ".*CPFEM.f90")
list (FILTER damask-sources EXCLUDE REGEX ".*DAMASK_marc.*.f90") list(FILTER damask-sources EXCLUDE REGEX ".*DAMASK_marc.*.f90")
list (FILTER damask-sources EXCLUDE REGEX ".*mesh_marc.*.f90") list(FILTER damask-sources EXCLUDE REGEX ".*mesh_marc.*.f90")
list (FILTER damask-sources EXCLUDE REGEX ".*mesh_abaqus.*.f90") list(FILTER damask-sources EXCLUDE REGEX ".*commercialFEM_fileList.*.f90")
list (FILTER damask-sources EXCLUDE REGEX ".*commercialFEM_fileList.*.f90")
if (PROJECT_NAME STREQUAL "damask-grid") if (PROJECT_NAME STREQUAL "damask-grid")
list (FILTER damask-sources EXCLUDE REGEX ".*mesh_FEM.*.f90") list(FILTER damask-sources EXCLUDE REGEX ".*mesh_FEM.*.f90")
file (GLOB grid-sources grid/*.f90) file(GLOB grid-sources grid/*.f90)
if (NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
add_executable (DAMASK_spectral ${damask-sources} ${grid-sources}) add_executable(DAMASK_spectral ${damask-sources} ${grid-sources})
install (TARGETS DAMASK_spectral RUNTIME DESTINATION bin) install (TARGETS DAMASK_spectral RUNTIME DESTINATION bin)
else() else()
add_library (DAMASK_spectral OBJECT ${damask-sources} ${grid-sources}) add_library(DAMASK_spectral OBJECT ${damask-sources} ${grid-sources})
exec_program (mktemp OUTPUT_VARIABLE nothing) exec_program (mktemp OUTPUT_VARIABLE nothing)
exec_program (mktemp ARGS -d OUTPUT_VARIABLE black_hole) exec_program (mktemp ARGS -d OUTPUT_VARIABLE black_hole)
install (PROGRAMS ${nothing} DESTINATION ${black_hole}) install (PROGRAMS ${nothing} DESTINATION ${black_hole})
endif() endif()
elseif (PROJECT_NAME STREQUAL "damask-mesh") elseif (PROJECT_NAME STREQUAL "damask-mesh")
list (FILTER damask-sources EXCLUDE REGEX ".*mesh_grid.*.f90") list(FILTER damask-sources EXCLUDE REGEX ".*mesh_grid.*.f90")
file (GLOB mesh-sources mesh/*.f90) file(GLOB mesh-sources mesh/*.f90)
add_executable (DAMASK_FEM ${damask-sources} ${mesh-sources}) add_executable(DAMASK_FEM ${damask-sources} ${mesh-sources})
install (TARGETS DAMASK_FEM RUNTIME DESTINATION bin) install (TARGETS DAMASK_FEM RUNTIME DESTINATION bin)
endif() endif()

3
src/CPFEM.f90
View File

@ -39,8 +39,7 @@ module CPFEM
integer(pInt), public :: & integer(pInt), public :: &
cycleCounter = 0_pInt, & !< needs description cycleCounter = 0_pInt, & !< needs description
theInc = -1_pInt, & !< needs description theInc = -1_pInt, & !< needs description
lastLovl = 0_pInt, & !< lovl in previous call to marc hypela2 lastLovl = 0_pInt !< lovl in previous call to marc hypela2
lastStep = 0_pInt !< kstep in previous call to abaqus umat
real(pReal), public :: & real(pReal), public :: &
theTime = 0.0_pReal, & !< needs description theTime = 0.0_pReal, & !< needs description
theDelta = 0.0_pReal theDelta = 0.0_pReal

380
src/DAMASK_abaqus.f
View File

@ -1,380 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Koen Janssens, Paul Scherrer Institut
!> @author Arun Prakash, Fraunhofer IWM
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief interfaces DAMASK with Abaqus/Standard
!> @details put the included file abaqus_v6.env in either your home or model directory,
!> it is a minimum Abaqus environment file containing all changes necessary to use the
!> DAMASK subroutine (see Abaqus documentation for more information on the use of abaqus_v6.env)
!> @details Abaqus subroutines used:
!> @details - UMAT
!> @details - DFLUX
!--------------------------------------------------------------------------------------------------
#define Abaqus
#include "prec.f90"
module DAMASK_interface
implicit none
private
character(len=4), dimension(2), parameter, public :: INPUTFILEEXTENSION = ['.pes','.inp']
public :: &
DAMASK_interface_init, &
getSolverJobName
contains
!--------------------------------------------------------------------------------------------------
!> @brief reports and sets working directory
!--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init
#if __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use ifport, only: &
CHDIR
implicit none
integer, dimension(8) :: &
dateAndTime
integer :: lenOutDir,ierr
character(len=256) :: wd
write(6,'(/,a)') ' <<<+- DAMASK_abaqus init -+>>>'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
write(6,'(/,a)') ' Version: '//DAMASKVERSION
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
#if __INTEL_COMPILER >= 1800
write(6,'(/,a)') ' Compiled with: '//compiler_version()
write(6,'(a)') ' Compiler options: '//compiler_options()
#else
write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
', build date :', __INTEL_COMPILER_BUILD_DATE
#endif
write(6,'(/,a)') ' Compiled on: '//__DATE__//' at '//__TIME__
call date_and_time(values = dateAndTime)
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)
call getoutdir(wd, lenOutDir)
ierr = CHDIR(wd)
if (ierr /= 0) then
write(6,'(a20,a,a16)') ' working directory "',trim(wd),'" does not exist'
call quit(1)
endif
end subroutine DAMASK_interface_init
!--------------------------------------------------------------------------------------------------
!> @brief using Abaqus/Standard function to get solver job name
!--------------------------------------------------------------------------------------------------
character(1024) function getSolverJobName()
implicit none
integer :: lenJobName
getSolverJobName=''
call getJobName(getSolverJobName, lenJobName)
end function getSolverJobName
end module DAMASK_interface
#include "commercialFEM_fileList.f90"
!--------------------------------------------------------------------------------------------------
!> @brief This is the Abaqus std user subroutine for defining material behavior
!--------------------------------------------------------------------------------------------------
subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,&
TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,NTENS,&
NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,&
DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
use prec, only: &
pReal, &
pInt
use numerics, only: &
!$ DAMASK_NumThreadsInt, &
usePingPong
use FEsolving, only: &
calcMode, &
terminallyIll, &
symmetricSolver
use debug, only: &
debug_info, &
debug_reset, &
debug_levelBasic, &
debug_level, &
debug_abaqus
use mesh, only: &
mesh_unitlength, &
mesh_FEasCP, &
mesh_ipCoordinates
use CPFEM, only: &
CPFEM_general, &
CPFEM_init_done, &
CPFEM_initAll, &
CPFEM_CALCRESULTS, &
CPFEM_AGERESULTS, &
CPFEM_COLLECT, &
CPFEM_BACKUPJACOBIAN, &
cycleCounter, &
theInc, &
theTime, &
theDelta, &
lastIncConverged, &
outdatedByNewInc, &
outdatedFFN1, &
lastStep
implicit none
integer(pInt), intent(in) :: &
nDi, & !< Number of direct stress components at this point
nShr, & !< Number of engineering shear stress components at this point
nTens, & !< Size of the stress or strain component array (NDI + NSHR)
nStatV, & !< Number of solution-dependent state variables
nProps, & !< User-defined number of material constants
noEl, & !< element number
nPt,& !< integration point number
kSlay, & !< layer number (shell elements etc.)
kSpt, & !< section point within the current layer
kStep, & !< step number
kInc !< increment number
character(len=80), intent(in) :: &
cmname !< uses-specified material name, left justified
real(pReal), intent(in) :: &
DTIME, &
TEMP, &
DTEMP, &
CELENT
real(pReal), dimension(1), intent(in) :: &
PREDEF, &
DPRED
real(pReal), dimension(2), intent(in) :: &
TIME !< step time/total time at beginning of the current increment
real(pReal), dimension(3), intent(in) :: &
COORDS
real(pReal), dimension(nTens), intent(in) :: &
STRAN, & !< total strains at beginning of the increment
DSTRAN !< strain increments
real(pReal), dimension(nProps), intent(in) :: &
PROPS
real(pReal), dimension(3,3), intent(in) :: &
DROT, & !< rotation increment matrix
DFGRD0, & !< F at beginning of increment
DFGRD1 !< F at end of increment
real(pReal), intent(inout) :: &
PNEWDT, & !< ratio of suggested new time increment
SSE, & !< specific elastic strain engergy
SPD, & !< specific plastic dissipation
SCD, & !< specific creep dissipation
RPL, & !< volumetric heat generation per unit time at the end of the increment
DRPLDT !< varation of RPL with respect to the temperature
real(pReal), dimension(nTens), intent(inout) :: &
STRESS !< stress tensor at the beginning of the increment, needs to be updated
real(pReal), dimension(nStatV), intent(inout) :: &
STATEV !< solution-dependent state variables
real(pReal), dimension(nTens), intent(out) :: &
DDSDDT, &
DRPLDE
real(pReal), dimension(nTens,nTens), intent(out) :: &
DDSDDE !< Jacobian matrix of the constitutive model
real(pReal) :: temperature ! temp by Abaqus is intent(in)
real(pReal), dimension(6) :: stress_h
real(pReal), dimension(6,6) :: ddsdde_h
integer(pInt) :: computationMode, i, cp_en
logical :: cutBack
#ifdef _OPENMP
integer :: defaultNumThreadsInt !< default value set by Abaqus
include "omp_lib.h"
defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
#endif
temperature = temp ! temp is intent(in)
DDSDDT = 0.0_pReal
DRPLDE = 0.0_pReal
if (iand(debug_level(debug_abaqus),debug_levelBasic) /= 0 .and. noel == 1 .and. npt == 1) then
write(6,*) 'el',noel,'ip',npt
write(6,*) 'got kInc as',kInc
write(6,*) 'got dStran',dStran
flush(6)
endif
if (.not. CPFEM_init_done) call CPFEM_initAll(noel,npt)
computationMode = 0
cp_en = mesh_FEasCP('elem',noel)
if (time(2) > theTime .or. kInc /= theInc) then ! reached convergence
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
if (kInc == 1) then ! >> start of analysis <<
lastIncConverged = .false. ! no Jacobian backup
outdatedByNewInc = .false. ! no aging of state
calcMode = .false. ! pretend last step was collection
write (6,'(i8,1x,i2,1x,a)') noel,npt,'<< UMAT >> start of analysis..!';flush(6)
else if (kInc - theInc > 1) then ! >> restart of broken analysis <<
lastIncConverged = .false. ! no Jacobian backup
outdatedByNewInc = .false. ! no aging of state
calcMode = .true. ! pretend last step was calculation
write (6,'(i8,1x,i2,1x,a)') noel,npt,'<< UMAT >> restart of analysis..!';flush(6)
else ! >> just the next inc <<
lastIncConverged = .true. ! request Jacobian backup
outdatedByNewInc = .true. ! request aging of state
calcMode = .true. ! assure last step was calculation
write (6,'(i8,1x,i2,1x,a)') noel,npt,'<< UMAT >> new increment..!';flush(6)
endif
else if ( dtime < theDelta ) then ! >> cutBack <<
lastIncConverged = .false. ! no Jacobian backup
outdatedByNewInc = .false. ! no aging of state
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
calcMode = .true. ! pretend last step was calculation
write(6,'(i8,1x,i2,1x,a)') noel,npt,'<< UMAT >> cutback detected..!';flush(6)
endif ! convergence treatment end
if (usePingPong) then
calcMode(npt,cp_en) = .not. calcMode(npt,cp_en) ! ping pong (calc <--> collect)
if (calcMode(npt,cp_en)) then ! now --- CALC ---
computationMode = CPFEM_CALCRESULTS
if ( lastStep /= kStep ) then ! first after ping pong
call debug_reset() ! resets debugging
outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1_pInt
endif
if(outdatedByNewInc) then
computationMode = ior(computationMode,CPFEM_AGERESULTS) ! calc and age results
outdatedByNewInc = .false. ! reset flag
endif
else ! now --- COLLECT ---
computationMode = CPFEM_COLLECT ! plain collect
if(lastStep /= kStep .and. .not. terminallyIll) &
call debug_info() ! first after ping pong reports (meaningful) debugging
if (lastIncConverged) then
computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! collect and backup Jacobian after convergence
lastIncConverged = .false. ! reset flag
endif
mesh_ipCoordinates(1:3,npt,cp_en) = mesh_unitlength * COORDS
endif
else ! --- PLAIN MODE ---
computationMode = CPFEM_CALCRESULTS ! always calc
if (lastStep /= kStep) then
if (.not. terminallyIll) &
call debug_info() ! first reports (meaningful) debugging
call debug_reset() ! and resets debugging
outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1_pInt
endif
if (outdatedByNewInc) then
computationMode = ior(computationMode,CPFEM_AGERESULTS)
outdatedByNewInc = .false. ! reset flag
endif
if (lastIncConverged) then
computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! backup Jacobian after convergence
lastIncConverged = .false. ! reset flag
endif
endif
theTime = time(2) ! record current starting time
theDelta = dtime ! record current time increment
theInc = kInc ! record current increment number
lastStep = kStep ! record step number
if (iand(debug_level(debug_abaqus),debug_levelBasic) /= 0) then
write(6,'(a16,1x,i2,1x,a,i8,a,i8,1x,i5,a)') 'computationMode',computationMode,'(',cp_en,':',noel,npt,')'
flush(6)
endif
call CPFEM_general(computationMode,usePingPong,dfgrd0,dfgrd1,temperature,dtime,noel,npt,stress_h,ddsdde_h)
! DAMASK: 11, 22, 33, 12, 23, 13
! ABAQUS explicit: 11, 22, 33, 12, 23, 13
! ABAQUS implicit: 11, 22, 33, 12, 13, 23
! ABAQUS implicit: 11, 22, 33, 12
ddsdde = ddsdde_h(1:ntens,1:ntens)
stress = stress_h(1:ntens)
if(symmetricSolver) ddsdde = 0.5_pReal*(ddsdde + transpose(ddsdde))
if(ntens == 6) then
stress_h = stress
stress(5) = stress_h(6)
stress(6) = stress_h(5)
ddsdde_h = ddsdde
ddsdde(:,5) = ddsdde_h(:,6)
ddsdde(:,6) = ddsdde_h(:,5)
ddsdde_h = ddsdde
ddsdde(5,:) = ddsdde_h(6,:)
ddsdde(6,:) = ddsdde_h(5,:)
end if
statev = 0
if (terminallyIll) pnewdt = 0.5_pReal ! force cutback directly ?
!$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value
end subroutine UMAT
!--------------------------------------------------------------------------------------------------
!> @brief calculate internal heat generated due to inelastic energy dissipation
!--------------------------------------------------------------------------------------------------
SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,&
JLTYP,TEMP,PRESS,SNAME)
use prec, only: &
pReal, &
pInt
use thermal_conduction, only: &
thermal_conduction_getSourceAndItsTangent
use mesh, only: &
mesh_FEasCP
implicit none
character(len=1024) :: sname
real(pReal), dimension(2), intent(out) :: flux
real(pReal), intent(in) :: sol
Integer(pInt), intent(in) :: Kstep, Kinc, Noel, Npt
real(pReal), dimension(3), intent(in) :: coords
real(pReal), intent(in) :: time
real(pReal) :: Jltyp
real(pReal), intent(in) :: temp
real(pReal), intent(in) :: Press
jltyp = 1
call thermal_conduction_getSourceAndItsTangent(flux(1), flux(2), sol, npt,mesh_FEasCP('elem',noel))
end subroutine DFLUX
!--------------------------------------------------------------------------------------------------
!> @brief calls the exit function of Abaqus/Standard
!--------------------------------------------------------------------------------------------------
subroutine quit(DAMASK_error)
use prec, only: &
pInt
implicit none
integer(pInt) :: DAMASK_error
flush(6)
call xit
end subroutine quit

2
src/FEsolving.f90
View File

@ -15,7 +15,7 @@ module FEsolving
FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element
FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP
#if defined(Marc4DAMASK) || defined(Abaqus) #if defined(Marc4DAMASK)
logical, dimension(:,:), allocatable :: & logical, dimension(:,:), allocatable :: &
calcMode !< do calculation or simply collect when using ping pong scheme calcMode !< do calculation or simply collect when using ping pong scheme
#endif #endif

242
src/IO.f90
View File

@ -33,16 +33,6 @@ module IO
IO_lc, & IO_lc, &
IO_error, & IO_error, &
IO_warning IO_warning
#if defined(Marc4DAMASK) || defined(Abaqus)
public :: &
IO_open_inputFile
#if defined(Abaqus)
public :: &
IO_continuousIntValues, &
IO_extractValue, &
IO_countDataLines
#endif
#endif
contains contains
@ -175,94 +165,6 @@ integer function IO_open_binary(fileName,mode)
end function IO_open_binary end function IO_open_binary
#if defined(Marc4DAMASK) || defined(Abaqus)
!--------------------------------------------------------------------------------------------------
!> @brief opens FEM input file for reading located in current working directory to given unit
!--------------------------------------------------------------------------------------------------
subroutine IO_open_inputFile(fileUnit)
integer, intent(in) :: fileUnit !< file unit
integer :: myStat
character(len=1024) :: path
#if defined(Abaqus)
integer :: fileType
fileType = 1 ! assume .pes
path = trim(getSolverJobName())//inputFileExtension(fileType) ! attempt .pes, if it exists: it should be used
open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
if(myStat /= 0) then ! if .pes does not work / exist; use conventional extension, i.e.".inp"
fileType = 2
path = trim(getSolverJobName())//inputFileExtension(fileType)
open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
endif
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path)
path = trim(getSolverJobName())//inputFileExtension(fileType)//'_assembly'
open(fileUnit,iostat=myStat,file=path)
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path)
if (.not.abaqus_assembleInputFile(fileUnit,fileUnit+1)) call IO_error(103) ! strip comments and concatenate any "include"s
close(fileUnit+1)
contains
!--------------------------------------------------------------------------------------------------
!> @brief create a new input file for abaqus simulations by removing all comment lines and
!> including "include"s
!--------------------------------------------------------------------------------------------------
recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess)
integer, intent(in) :: unit1, &
unit2
integer, allocatable, dimension(:) :: chunkPos
character(len=pStringLen :: line,fname
logical :: createSuccess,fexist
do
read(unit2,'(A)',END=220) line
chunkPos = IO_stringPos(line)
if (IO_lc(IO_StringValue(line,chunkPos,1))=='*include') then
fname = trim(line(9+scan(line(9:),'='):))
inquire(file=fname, exist=fexist)
if (.not.(fexist)) then
write(6,*)'ERROR: file does not exist error in abaqus_assembleInputFile'
write(6,*)'filename: ', trim(fname)
createSuccess = .false.
return
endif
open(unit2+1,err=200,status='old',file=fname)
if (abaqus_assembleInputFile(unit1,unit2+1)) then
createSuccess=.true.
close(unit2+1)
else
createSuccess=.false.
return
endif
else if (line(1:2) /= '**' .OR. line(1:8)=='**damask') then
write(unit1,'(A)') trim(line)
endif
enddo
220 createSuccess = .true.
return
200 createSuccess =.false.
end function abaqus_assembleInputFile
#elif defined(Marc4DAMASK)
path = trim(getSolverJobName())//inputFileExtension
open(fileUnit,status='old',iostat=myStat,file=path)
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path)
#endif
end subroutine IO_open_inputFile
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief identifies strings without content !> @brief identifies strings without content
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -419,7 +321,6 @@ end function IO_lc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx !> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx
!> in ABAQUS either time step is reduced or execution terminated
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
@ -608,29 +509,6 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
case (894) case (894)
msg = 'MPI error' msg = 'MPI error'
!-------------------------------------------------------------------------------------------------
! error messages related to parsing of Abaqus input file
case (900)
msg = 'improper definition of nodes in input file (Nnodes < 2)'
case (901)
msg = 'no elements defined in input file (Nelems = 0)'
case (902)
msg = 'no element sets defined in input file (No *Elset exists)'
case (903)
msg = 'no materials defined in input file (Look into section assigments)'
case (904)
msg = 'no elements could be assigned for Elset: '
case (905)
msg = 'error in mesh_abaqus_map_materials'
case (906)
msg = 'error in mesh_abaqus_count_cpElements'
case (907)
msg = 'size of mesh_mapFEtoCPelem in mesh_abaqus_map_elements'
case (908)
msg = 'size of mesh_mapFEtoCPnode in mesh_abaqus_map_nodes'
case (909)
msg = 'size of mesh_node in mesh_abaqus_build_nodes not equal to mesh_Nnodes'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! general error messages ! general error messages
@ -751,126 +629,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
end subroutine IO_warning end subroutine IO_warning
#ifdef Abaqus
!--------------------------------------------------------------------------------------------------
!> @brief extracts string value from key=value pair and check whether key matches
!--------------------------------------------------------------------------------------------------
character(len=300) pure function IO_extractValue(pair,key)
character(len=*), intent(in) :: pair, & !< key=value pair
key !< key to be expected
character(len=*), parameter :: SEP = achar(61) ! '='
integer :: myChunk !< position number of desired chunk
IO_extractValue = ''
myChunk = scan(pair,SEP)
if (myChunk > 0 .and. pair(:myChunk-1) == key) IO_extractValue = pair(myChunk+1:) ! extract value if key matches
end function IO_extractValue
!--------------------------------------------------------------------------------------------------
!> @brief count lines containig data up to next *keyword
!--------------------------------------------------------------------------------------------------
integer function IO_countDataLines(fileUnit)
integer, intent(in) :: fileUnit !< file handle
integer, allocatable, dimension(:) :: chunkPos
character(len=pStringLen) :: line, &
tmp
IO_countDataLines = 0
line = ''
do while (trim(line) /= IO_EOF)
read(fileUnit,'(A)') line
chunkPos = IO_stringPos(line)
tmp = IO_lc(IO_stringValue(line,chunkPos,1))
if (tmp(1:1) == '*' .and. tmp(2:2) /= '*') then ! found keyword
exit
else
if (tmp(2:2) /= '*') IO_countDataLines = IO_countDataLines + 1
endif
enddo
backspace(fileUnit)
end function IO_countDataLines
!--------------------------------------------------------------------------------------------------
!> @brief return integer list corresponding to items in consecutive lines.
!! First integer in array is counter
!> @details Abaqus: triplet of start,stop,inc or named set
!--------------------------------------------------------------------------------------------------
function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
integer, intent(in) :: maxN
integer, dimension(1+maxN) :: IO_continuousIntValues
integer, intent(in) :: fileUnit, &
lookupMaxN
integer, dimension(:,:), intent(in) :: lookupMap
character(len=*), dimension(:), intent(in) :: lookupName
integer :: i,first,last
integer :: j,l,c
integer, allocatable, dimension(:) :: chunkPos
character(len=pStringLen) :: line
logical :: rangeGeneration
IO_continuousIntValues = 0
rangeGeneration = .false.
c = IO_countDataLines(fileUnit)
do l = 1,c
backspace(fileUnit)
enddo
!--------------------------------------------------------------------------------------------------
! check if the element values in the elset are auto generated
backspace(fileUnit)
read(fileUnit,'(A)',end=100) line
chunkPos = IO_stringPos(line)
do i = 1,chunkPos(1)
if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'generate') rangeGeneration = .true.
enddo
do l = 1,c
read(fileUnit,'(A)',end=100) line
chunkPos = IO_stringPos(line)
if (verify(IO_stringValue(line,chunkPos,1),'0123456789') > 0) then ! a non-int, i.e. set names follow on this line
do i = 1,chunkPos(1) ! loop over set names in line
do j = 1,lookupMaxN ! look through known set names
if (IO_stringValue(line,chunkPos,i) == lookupName(j)) then ! found matching name
first = 2 + IO_continuousIntValues(1) ! where to start appending data
last = first + lookupMap(1,j) - 1 ! up to where to append data
IO_continuousIntValues(first:last) = lookupMap(2:1+lookupMap(1,j),j) ! add resp. entity list
IO_continuousIntValues(1) = IO_continuousIntValues(1) + lookupMap(1,j) ! count them
endif
enddo
enddo
else if (rangeGeneration) then ! range generation
do i = IO_intValue(line,chunkPos,1),&
IO_intValue(line,chunkPos,2),&
max(1,IO_intValue(line,chunkPos,3))
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1
IO_continuousIntValues(1+IO_continuousIntValues(1)) = i
enddo
else ! read individual elem nums
do i = 1,chunkPos(1)
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1
IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,chunkPos,i)
enddo
endif
enddo
100 end function IO_continuousIntValues
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! internal helper functions ! internal helper functions

5
src/commercialFEM_fileList.f90
View File

@ -1,7 +1,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief all DAMASK files without solver !> @brief all DAMASK files without solver
!> @details List of files needed by MSC.Marc and Abaqus/Standard !> @details List of files needed by MSC.Marc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
#include "IO.f90" #include "IO.f90"
#include "numerics.f90" #include "numerics.f90"
@ -19,9 +19,6 @@
#include "results.f90" #include "results.f90"
#include "geometry_plastic_nonlocal.f90" #include "geometry_plastic_nonlocal.f90"
#include "discretization.f90" #include "discretization.f90"
#ifdef Abaqus
#include "mesh_abaqus.f90"
#endif
#ifdef Marc4DAMASK #ifdef Marc4DAMASK
#include "mesh_marc.f90" #include "mesh_marc.f90"
#endif #endif

20
src/constitutive.f90
View File

@ -202,7 +202,7 @@ module constitutive
of of
end subroutine plastic_disloUCLA_dotState end subroutine plastic_disloUCLA_dotState
module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, &
timestep,ip,el) timestep,ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< current integration point ip, & !< current integration point
@ -213,7 +213,7 @@ module constitutive
real(pReal), dimension(3,3), intent(in) ::& real(pReal), dimension(3,3), intent(in) ::&
Mp !< MandelStress Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: &
Fe, & !< elastic deformation gradient F, & !< deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
end subroutine plastic_nonlocal_dotState end subroutine plastic_nonlocal_dotState
@ -232,12 +232,12 @@ module constitutive
of of
end subroutine plastic_disloUCLA_dependentState end subroutine plastic_disloUCLA_dependentState
module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) module subroutine plastic_nonlocal_dependentState(F, Fp, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & ip, &
el el
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Fe, & F, &
Fp Fp
end subroutine plastic_nonlocal_dependentState end subroutine plastic_nonlocal_dependentState
@ -412,14 +412,14 @@ end function constitutive_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different constitutive models !> @brief calls microstructure function of the different constitutive models
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) subroutine constitutive_microstructure(F, Fp, ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient F, & !< elastic deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
integer :: & integer :: &
ho, & !< homogenization ho, & !< homogenization
@ -439,7 +439,7 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el)) instance = phase_plasticityInstance(material_phaseAt(ipc,el))
call plastic_disloUCLA_dependentState(instance,of) call plastic_disloUCLA_dependentState(instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dependentState (Fe,Fp,ip,el) call plastic_nonlocal_dependentState (F,Fp,ip,el)
end select plasticityType end select plasticityType
end subroutine constitutive_microstructure end subroutine constitutive_microstructure
@ -715,7 +715,7 @@ end subroutine constitutive_hooke_SandItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el) subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
@ -724,7 +724,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
FeArray, & !< elastic deformation gradient FArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
@ -771,7 +771,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (Mp,FeArray,FpArray,temperature(ho)%p(tme), & call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme), &
subdt,ip,el) subdt,ip,el)
end select plasticityType end select plasticityType

125
src/constitutive_plastic_nonlocal.f90
View File

@ -180,7 +180,8 @@ submodule(constitutive) plastic_nonlocal
type(tNonlocalState), allocatable, dimension(:) :: & type(tNonlocalState), allocatable, dimension(:) :: &
deltaState, & deltaState, &
dotState, & dotState, &
state state, &
state0
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
@ -226,6 +227,7 @@ module subroutine plastic_nonlocal_init
allocate(param(maxNinstances)) allocate(param(maxNinstances))
allocate(state(maxNinstances)) allocate(state(maxNinstances))
allocate(state0(maxNinstances))
allocate(dotState(maxNinstances)) allocate(dotState(maxNinstances))
allocate(deltaState(maxNinstances)) allocate(deltaState(maxNinstances))
allocate(microstructure(maxNinstances)) allocate(microstructure(maxNinstances))
@ -238,6 +240,7 @@ module subroutine plastic_nonlocal_init
associate(prm => param(phase_plasticityInstance(p)), & associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), &
st0 => state0(phase_plasticityInstance(p)), &
del => deltaState(phase_plasticityInstance(p)), & del => deltaState(phase_plasticityInstance(p)), &
dst => microstructure(phase_plasticityInstance(p)), & dst => microstructure(phase_plasticityInstance(p)), &
config => config_phase(p)) config => config_phase(p))
@ -482,6 +485,7 @@ module subroutine plastic_nonlocal_init
totalNslip(phase_plasticityInstance(p)) = prm%totalNslip totalNslip(phase_plasticityInstance(p)) = prm%totalNslip
st0%rho => plasticState(p)%state0 (0*prm%totalNslip+1:10*prm%totalNslip,:)
stt%rho => plasticState(p)%state (0*prm%totalNslip+1:10*prm%totalNslip,:) stt%rho => plasticState(p)%state (0*prm%totalNslip+1:10*prm%totalNslip,:)
dot%rho => plasticState(p)%dotState (0*prm%totalNslip+1:10*prm%totalNslip,:) dot%rho => plasticState(p)%dotState (0*prm%totalNslip+1:10*prm%totalNslip,:)
del%rho => plasticState(p)%deltaState (0*prm%totalNslip+1:10*prm%totalNslip,:) del%rho => plasticState(p)%deltaState (0*prm%totalNslip+1:10*prm%totalNslip,:)
@ -712,13 +716,13 @@ end subroutine plastic_nonlocal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure !> @brief calculates quantities characterizing the microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) module subroutine plastic_nonlocal_dependentState(F, Fp, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & ip, &
el el
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Fe, & F, &
Fp Fp
integer :: & integer :: &
@ -761,7 +765,8 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
rho_scr_delta rho_scr_delta
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
rho, & rho, &
rho_neighbor rho0, &
rho_neighbor0
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))), & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))), &
totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
myInteractionMatrix ! corrected slip interaction matrix myInteractionMatrix ! corrected slip interaction matrix
@ -791,7 +796,7 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
! coefficients are corrected for the line tension effect ! coefficients are corrected for the line tension effect
! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals)
if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then ! only fcc and bcc if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then
do s = 1,ns do s = 1,ns
correction = ( 1.0_pReal - prm%linetensionEffect & correction = ( 1.0_pReal - prm%linetensionEffect &
+ prm%linetensionEffect & + prm%linetensionEffect &
@ -815,13 +820,13 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
! ToDo: MD: this is most likely only correct for F_i = I ! ToDo: MD: this is most likely only correct for F_i = I
!################################################################################################# !#################################################################################################
rho0 = getRho0(instance,of,ip,el)
if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then
invFe = math_inv33(Fe)
invFp = math_inv33(Fp) invFp = math_inv33(Fp)
invFe = matmul(Fp,math_inv33(F))
rho_edg_delta = rho(:,mob_edg_pos) - rho(:,mob_edg_neg) rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
rho_scr_delta = rho(:,mob_scr_pos) - rho(:,mob_scr_neg) rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
rhoExcess(1,1:ns) = rho_edg_delta rhoExcess(1,1:ns) = rho_edg_delta
rhoExcess(2,1:ns) = rho_scr_delta rhoExcess(2,1:ns) = rho_scr_delta
@ -841,13 +846,13 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
if (neighbor_instance == instance) then if (neighbor_instance == instance) then
nRealNeighbors = nRealNeighbors + 1.0_pReal nRealNeighbors = nRealNeighbors + 1.0_pReal
rho_neighbor = getRho(instance,no,neighbor_ip,neighbor_el) rho_neighbor0 = getRho0(instance,no,neighbor_ip,neighbor_el)
rho_edg_delta_neighbor(:,n) = rho_neighbor(:,mob_edg_pos) - rho_neighbor(:,mob_edg_neg) rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg)
rho_scr_delta_neighbor(:,n) = rho_neighbor(:,mob_scr_pos) - rho_neighbor(:,mob_scr_neg) rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg)
neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor(:,edg)),2) neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2)
neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor(:,scr)),2) neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2)
connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) & connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) &
- discretization_IPcoords(1:3,el+neighbor_ip-1)) - discretization_IPcoords(1:3,el+neighbor_ip-1))
@ -1320,7 +1325,7 @@ end subroutine plastic_nonlocal_deltaState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, &
timestep,ip,el) timestep,ip,el)
integer, intent(in) :: & integer, intent(in) :: &
@ -1332,7 +1337,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
real(pReal), dimension(3,3), intent(in) ::& real(pReal), dimension(3,3), intent(in) ::&
Mp !< MandelStress Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: &
Fe, & !< elastic deformation gradient F, & !< elastic deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
integer :: & integer :: &
@ -1358,6 +1363,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
s !< index of my current slip system s !< index of my current slip system
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
rho, & rho, &
rho0, & !< dislocation density at beginning of time step
rhoDot, & !< density evolution rhoDot, & !< density evolution
rhoDotMultiplication, & !< density evolution by multiplication rhoDotMultiplication, & !< density evolution by multiplication
rhoDotFlux, & !< density evolution by flux rhoDotFlux, & !< density evolution by flux
@ -1366,12 +1372,12 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
rhoDotThermalAnnihilation !< density evolution by thermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
my_v, & !< dislocation glide velocity of central ip v0, &
neighbor_v, & !< dislocation glide velocity of enighboring ip neighbor_v0, & !< dislocation glide velocity of enighboring ip
gdot !< shear rates gdot !< shear rates
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
tau, & !< current resolved shear stress tau, & !< current resolved shear stress
@ -1423,9 +1429,11 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
tau = 0.0_pReal tau = 0.0_pReal
gdot = 0.0_pReal gdot = 0.0_pReal
rho = getRho(instance,o,ip,el) rho = getRho(instance,o,ip,el)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
rhoDip = rho(:,dip) rhoDip = rho(:,dip)
rho0 = getRho0(instance,o,ip,el)
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) forall (s = 1:ns, t = 1:4)
v(s,t) = plasticState(p)%state(iV (s,t,instance),o) v(s,t) = plasticState(p)%state(iV (s,t,instance),o)
@ -1488,6 +1496,9 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
* sqrt(stt%rho_forest(:,o)) / prm%lambda0 / prm%burgers(1:ns), 2, 4) * sqrt(stt%rho_forest(:,o)) / prm%lambda0 / prm%burgers(1:ns), 2, 4)
endif isBCC endif isBCC
forall (s = 1:ns, t = 1:4)
v0(s,t) = plasticState(p)%state0(iV(s,t,instance),o)
endforall
!**************************************************************************** !****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal plasticity) !*** calculate dislocation fluxes (only for nonlocal plasticity)
@ -1496,14 +1507,14 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... if (any( abs(gdot) > 0.0_pReal & ! any active slip system ...
.and. prm%CFLfactor * abs(v) * timestep & .and. prm%CFLfactor * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then
write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', &
maxval(abs(v), abs(gdot) > 0.0_pReal & maxval(abs(v0), abs(gdot) > 0.0_pReal &
.and. prm%CFLfactor * abs(v) * timestep & .and. prm%CFLfactor * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
' at a timestep of ',timestep ' at a timestep of ',timestep
write(6,'(a)') '<< CONST >> enforcing cutback !!!' write(6,'(a)') '<< CONST >> enforcing cutback !!!'
@ -1522,8 +1533,8 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
m(1:3,1:ns,3) = -prm%slip_transverse m(1:3,1:ns,3) = -prm%slip_transverse
m(1:3,1:ns,4) = prm%slip_transverse m(1:3,1:ns,4) = prm%slip_transverse
my_Fe = Fe(1:3,1:3,1,ip,el) my_F = F(1:3,1:3,1,ip,el)
my_F = matmul(my_Fe, Fp(1:3,1:3,1,ip,el)) my_Fe = matmul(my_F, math_inv33(Fp(1:3,1:3,1,ip,el)))
neighbors: do n = 1,nIPneighbors neighbors: do n = 1,nIPneighbors
@ -1540,8 +1551,8 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el) neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el)
neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el)) neighbor_Fe = matmul(neighbor_F, math_inv33(Fp(1:3,1:3,1,neighbor_ip,neighbor_el)))
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average
Favg = my_F Favg = my_F
@ -1558,8 +1569,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
!* compatibility !* compatibility
considerEnteringFlux = .false. considerEnteringFlux = .false.
neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below
neighbor_rhoSgl = 0.0_pReal
if (neighbor_n > 0) then if (neighbor_n > 0) then
if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID & if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID &
.and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) & .and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) &
@ -1568,14 +1578,14 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
enteringFlux: if (considerEnteringFlux) then enteringFlux: if (considerEnteringFlux) then
forall (s = 1:ns, t = 1:4) forall (s = 1:ns, t = 1:4)
neighbor_v(s,t) = plasticState(np)%state(iV (s,t,neighbor_instance),no) neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no)
neighbor_rhoSgl(s,t) = max(plasticState(np)%state(iRhoU(s,t,neighbor_instance),no), & neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no), &
0.0_pReal) 0.0_pReal)
endforall endforall
where (neighbor_rhoSgl * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN &
.or. neighbor_rhoSgl < prm%significantRho) & .or. neighbor_rhoSgl0 < prm%significantRho) &
neighbor_rhoSgl = 0.0_pReal neighbor_rhoSgl0 = 0.0_pReal
normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!) IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) & normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
@ -1586,9 +1596,9 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
do t = 1,4 do t = 1,4
c = (t + 1) / 2 c = (t + 1) / 2
topp = t + mod(t,2) - mod(t+1,2) topp = t + mod(t,2) - mod(t+1,2)
if (neighbor_v(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me
.and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density
lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) &
* math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility...
rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) & rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) &
@ -1619,9 +1629,6 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
endif endif
leavingFlux: if (considerLeavingFlux) then leavingFlux: if (considerLeavingFlux) then
my_rhoSgl = rhoSgl
my_v = v
normal_me2neighbor_defConf = math_det33(Favg) & normal_me2neighbor_defConf = math_det33(Favg) &
* matmul(math_inv33(transpose(Favg)), & * matmul(math_inv33(transpose(Favg)), &
IPareaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) IPareaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!)
@ -1632,18 +1639,18 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
do s = 1,ns do s = 1,ns
do t = 1,4 do t = 1,4
c = (t + 1) / 2 c = (t + 1) / 2
if (my_v(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
if (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then ! no sign change in flux density if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density
transmissivity = sum(compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor transmissivity = sum(compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal transmissivity = 0.0_pReal
endif endif
lineLength = my_rhoSgl(s,t) * my_v(s,t) & lineLength = my_rhoSgl0(s,t) * v0(s,t) &
* math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
+ lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) &
* sign(1.0_pReal, my_v(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point * sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
endif endif
enddo enddo
enddo enddo
@ -1713,7 +1720,6 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = 0.0_pReal
rhoDot = rhoDotFlux & rhoDot = rhoDotFlux &
+ rhoDotMultiplication & + rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide & + rhoDotSingle2DipoleGlide &
@ -1931,6 +1937,31 @@ function getRho(instance,of,ip,el)
end function getRho end function getRho
!--------------------------------------------------------------------------------------------------
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
function getRho0(instance,of,ip,el)
integer, intent(in) :: instance, of,ip,el
real(pReal), dimension(param(instance)%totalNslip,10) :: getRho0
associate(prm => param(instance))
getRho0 = reshape(state0(instance)%rho(:,of),[prm%totalNslip,10])
! ensure positive densities (not for imm, they have a sign)
getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)
getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal)
where(abs(getRho0) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) &
getRho0 = 0.0_pReal
end associate
end function getRho0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

8
src/crystallite.f90
View File

@ -264,8 +264,8 @@ subroutine crystallite_init
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), & call constitutive_microstructure(crystallite_partionedF0(1:3,1:3,c,i,e), &
crystallite_Fp(1:3,1:3,c,i,e), & crystallite_partionedFp0(1:3,1:3,c,i,e), &
c,i,e) ! update dependent state variables to be consistent with basic states c,i,e) ! update dependent state variables to be consistent with basic states
enddo enddo
enddo enddo
@ -1963,9 +1963,9 @@ subroutine update_dotState(timeFraction)
!$OMP FLUSH(nonlocalStop) !$OMP FLUSH(nonlocalStop)
if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_Fe, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e)*timeFraction, g,i,e) crystallite_subdt(g,i,e)*timeFraction, g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))

9
src/debug.f90
View File

@ -37,10 +37,9 @@ module debug
debug_HOMOGENIZATION = 9, & debug_HOMOGENIZATION = 9, &
debug_CPFEM = 10, & debug_CPFEM = 10, &
debug_SPECTRAL = 11, & debug_SPECTRAL = 11, &
debug_MARC = 12, & debug_MARC = 12
debug_ABAQUS = 13
integer, parameter, private :: & integer, parameter, private :: &
debug_MAXNTYPE = debug_ABAQUS !< must be set to the maximum defined debug type debug_MAXNTYPE = debug_MARC !< must be set to the maximum defined debug type
integer,protected, dimension(debug_maxNtype+2), public :: & ! specific ones, and 2 for "all" and "other" integer,protected, dimension(debug_maxNtype+2), public :: & ! specific ones, and 2 for "all" and "other"
debug_level = 0 debug_level = 0
@ -136,8 +135,6 @@ subroutine debug_init
what = debug_SPECTRAL what = debug_SPECTRAL
case ('marc') case ('marc')
what = debug_MARC what = debug_MARC
case ('abaqus')
what = debug_ABAQUS
case ('all') case ('all')
what = debug_MAXNTYPE + 1 what = debug_MAXNTYPE + 1
case ('other') case ('other')
@ -210,8 +207,6 @@ subroutine debug_init
tag = ' Spectral solver' tag = ' Spectral solver'
case (debug_MARC) case (debug_MARC)
tag = ' MSC.MARC FEM solver' tag = ' MSC.MARC FEM solver'
case (debug_ABAQUS)
tag = ' ABAQUS FEM solver'
end select end select
if(debug_level(i) /= 0) then if(debug_level(i) /= 0) then

2822
src/mesh_abaqus.f90
View File

File diff suppressed because it is too large Load Diff

4
src/prec.f90
View File

@ -12,11 +12,7 @@ module prec
implicit none implicit none
public public
! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds
#ifdef Abaqus
integer, parameter, public :: pReal = selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit)
#else
integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit) integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit)
#endif
#if(INT==8) #if(INT==8)
integer, parameter, public :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) integer, parameter, public :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit)
#else #else