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

Merge branch 'development' into 25_materialConfigPythonInclude

This commit is contained in:
Vitesh Shah 2018-06-06 15:34:08 +02:00
parent e1c5168040 1ff02d9712
commit 73e42a4479
13 changed files with 231 additions and 337 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
CMakeLists.txt
View File

@ -5,8 +5,8 @@ cmake_minimum_required (VERSION 2.8.8 FATAL_ERROR)
#--------------------------------------------------------------------------------------- #---------------------------------------------------------------------------------------
# Find PETSc from system environment # Find PETSc from system environment
set(PETSC_DIR $ENV{PETSC_DIR}) set(PETSC_DIR $ENV{PETSC_DIR})
if ("${PETSC_DIR}" STREQUAL "") if (PETSC_DIR STREQUAL "")
message (FATAL_ERROR "PETSC_DIR is not defined") message (FATAL_ERROR "PETSc location (PETSC_DIR) is not defined")
endif () endif ()
set (petsc_conf_variables "${PETSC_DIR}/lib/petsc/conf/variables") set (petsc_conf_variables "${PETSC_DIR}/lib/petsc/conf/variables")
@ -105,52 +105,54 @@ set (CMAKE_C_COMPILER "${PETSC_MPICC}")
# Now start to care about DAMASK # Now start to care about DAMASK
# DAMASK solver defines project to build # DAMASK solver defines project to build
if ("${DAMASK_SOLVER}" STREQUAL "SPECTRAL") if (DAMASK_SOLVER STREQUAL "SPECTRAL")
project (DAMASK_spectral Fortran C) project (DAMASK_spectral Fortran C)
add_definitions (-DSpectral) add_definitions (-DSpectral)
message ("Building Spectral Solver\n") message ("Building Spectral Solver\n")
elseif ("${DAMASK_SOLVER}" STREQUAL "FEM") elseif (DAMASK_SOLVER STREQUAL "FEM")
project (DAMASK_FEM Fortran C) project (DAMASK_FEM Fortran C)
add_definitions (-DFEM) add_definitions (-DFEM)
message ("Building FEM Solver\n") message ("Building FEM Solver\n")
else ()
message (FATAL_ERROR "Build target (DAMASK_SOLVER) is not defined")
endif () endif ()
# set linker commands (needs to be done after defining the project) # set linker commands (needs to be done after defining the project)
set (CMAKE_LINKER "${PETSC_LINKER}") set (CMAKE_LINKER "${PETSC_LINKER}")
if ("${CMAKE_BUILD_TYPE}" STREQUAL "") if (CMAKE_BUILD_TYPE STREQUAL "")
set (CMAKE_BUILD_TYPE "RELEASE") set (CMAKE_BUILD_TYPE "RELEASE")
endif () endif ()
# Predefined sets for OPTIMIZATION/OPENMP based on BUILD_TYPE # Predefined sets for OPTIMIZATION/OPENMP based on BUILD_TYPE
if ("${CMAKE_BUILD_TYPE}" STREQUAL "DEBUG" OR "${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY" ) if (CMAKE_BUILD_TYPE STREQUAL "DEBUG" OR CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -DDEBUG") set (DEBUG_FLAGS "${DEBUG_FLAGS} -DDEBUG")
set (PARALLEL "OFF") set (PARALLEL "OFF")
set (OPTI "OFF") set (OPTI "OFF")
elseif ("${CMAKE_BUILD_TYPE}" STREQUAL "RELEASE") elseif (CMAKE_BUILD_TYPE STREQUAL "RELEASE")
set (PARALLEL "ON") set (PARALLEL "ON")
set (OPTI "DEFENSIVE") set (OPTI "DEFENSIVE")
elseif ("${CMAKE_BUILD_TYPE}" STREQUAL "PERFORMANCE") elseif (CMAKE_BUILD_TYPE STREQUAL "PERFORMANCE")
set (PARALLEL "ON") set (PARALLEL "ON")
set (OPTI "AGGRESSIVE") set (OPTI "AGGRESSIVE")
endif () endif ()
# $OPTIMIZATION takes precedence over $BUILD_TYPE defaults # $OPTIMIZATION takes precedence over $BUILD_TYPE defaults
if ("${OPTIMIZATION}" STREQUAL "") if (OPTIMIZATION STREQUAL "")
set (OPTIMIZATION "${OPTI}") set (OPTIMIZATION "${OPTI}")
else () else ()
set (OPTIMIZATION "${OPTIMIZATION}") set (OPTIMIZATION "${OPTIMIZATION}")
endif () endif ()
# $OPENMP takes precedence over $BUILD_TYPE defaults # $OPENMP takes precedence over $BUILD_TYPE defaults
if ("${OPENMP}" STREQUAL "") if (OPENMP STREQUAL "")
set (OPENMP "${PARALLEL}") set (OPENMP "${PARALLEL}")
else () else ()
set(OPENMP "${OPENMP}") set(OPENMP "${OPENMP}")
endif () endif ()
# syntax check only (mainly for pre-receive hook, works only with gfortran) # syntax check only (mainly for pre-receive hook, works only with gfortran)
if ("${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY" ) if (CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
set (BUILDCMD_POST "${BUILDCMD_POST} -fsyntax-only") set (BUILDCMD_POST "${BUILDCMD_POST} -fsyntax-only")
endif () endif ()
@ -188,19 +190,19 @@ set (DAMASK_INCLUDE_FLAGS "${DAMASK_INCLUDE_FLAGS} ${PETSC_INCLUDES}")
################################################################################################### ###################################################################################################
# Intel Compiler # Intel Compiler
################################################################################################### ###################################################################################################
if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel") if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
if (OPENMP) if (OPENMP)
set (OPENMP_FLAGS "-qopenmp -parallel") set (OPENMP_FLAGS "-qopenmp -parallel")
endif () endif ()
if ("${OPTIMIZATION}" STREQUAL "OFF") if (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0 -no-ip") set (OPTIMIZATION_FLAGS "-O0 -no-ip")
elseif ("${OPTIMIZATION}" STREQUAL "DEFENSIVE") elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2") set (OPTIMIZATION_FLAGS "-O2")
elseif ("${OPTIMIZATION}" STREQUAL "AGGRESSIVE") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost") set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost")
# -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost" # -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost"
endif () endif ()
# -assume std_mod_proc_name (included in -standard-semantics) causes problems if other modules # -assume std_mod_proc_name (included in -standard-semantics) causes problems if other modules
@ -308,18 +310,18 @@ if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel")
################################################################################################### ###################################################################################################
# GNU Compiler # GNU Compiler
################################################################################################### ###################################################################################################
elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU") elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
if (OPENMP) if (OPENMP)
set (OPENMP_FLAGS "-fopenmp") set (OPENMP_FLAGS "-fopenmp")
endif () endif ()
if ("${OPTIMIZATION}" STREQUAL "OFF") if (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0" ) set (OPTIMIZATION_FLAGS "-O0" )
elseif ("${OPTIMIZATION}" STREQUAL "DEFENSIVE") elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2") set (OPTIMIZATION_FLAGS "-O2")
elseif ("${OPTIMIZATION}" STREQUAL "AGGRESSIVE") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-O3 -ffast-math -funroll-loops -ftree-vectorize") set (OPTIMIZATION_FLAGS "-O3 -ffast-math -funroll-loops -ftree-vectorize")
endif () endif ()
set (STANDARD_CHECK "-std=f2008ts -pedantic-errors" ) set (STANDARD_CHECK "-std=f2008ts -pedantic-errors" )
@ -443,12 +445,15 @@ elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU")
# Additional options # Additional options
# -fdefault-integer-8: Use it to set precision to 8 bytes for integer, don't use it for the standard case of pInt=4 (there is no -fdefault-integer-4) # -fdefault-integer-8: Use it to set precision to 8 bytes for integer, don't use it for the standard case of pInt=4 (there is no -fdefault-integer-4)
else ()
message (FATAL_ERROR "Compiler type (CMAKE_Fortran_COMPILER_ID) not recognized")
endif () endif ()
set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}") set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}")
set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${CMAKE_LINKER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}") set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${CMAKE_LINKER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}")
if ("${CMAKE_BUILD_TYPE}" STREQUAL "DEBUG") if (CMAKE_BUILD_TYPE STREQUAL "DEBUG")
set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}") set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}")
set (CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}") set (CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}")
endif () endif ()
@ -464,15 +469,15 @@ message ("Fortran Linker Command:\n${CMAKE_Fortran_LINK_EXECUTABLE}\n")
add_subdirectory (src) add_subdirectory (src)
# INSTALL BUILT BINARIES # INSTALL BUILT BINARIES
if ("${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY") if (CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
exec_program (mktemp ARGS -d OUTPUT_VARIABLE BLACK_HOLE) exec_program (mktemp ARGS -d OUTPUT_VARIABLE BLACK_HOLE)
install (PROGRAMS ${PROJECT_BINARY_DIR}/src/prec.mod install (PROGRAMS ${PROJECT_BINARY_DIR}/src/prec.mod
DESTINATION ${BLACK_HOLE}) DESTINATION ${BLACK_HOLE})
else () else ()
if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral") if (PROJECT_NAME STREQUAL "DAMASK_spectral")
install (PROGRAMS ${PROJECT_BINARY_DIR}/src/DAMASK_spectral install (PROGRAMS ${PROJECT_BINARY_DIR}/src/DAMASK_spectral
DESTINATION ${CMAKE_INSTALL_PREFIX}) DESTINATION ${CMAKE_INSTALL_PREFIX})
elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM") elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
install (PROGRAMS ${PROJECT_BINARY_DIR}/src/DAMASK_FEM install (PROGRAMS ${PROJECT_BINARY_DIR}/src/DAMASK_FEM
DESTINATION ${CMAKE_INSTALL_PREFIX}) DESTINATION ${CMAKE_INSTALL_PREFIX})
endif () endif ()

1
Makefile
View File

@ -25,7 +25,6 @@ build/FEM:
.PHONY: marc .PHONY: marc
marc: marc:
@./installation/symLink_Code.sh
@./installation/mods_MarcMentat/apply_DAMASK_modifications.sh ${MAKEFLAGS} @./installation/mods_MarcMentat/apply_DAMASK_modifications.sh ${MAKEFLAGS}
.PHONY: clean .PHONY: clean

2
PRIVATE

@ -1 +1 @@
Subproject commit 8a2e3ccce34be2a19cdcaecb8456a5dc6c2564f0 Subproject commit cd02f6c1a481491eb4517651516b8311348b4777

2
VERSION
View File

@ -1 +1 @@
v2.0.2-4-g35ba260 v2.0.2-37-gc903880

18
src/CMakeLists.txt
View File

@ -1,5 +1,5 @@
# special flags for some files # special flags for some files
if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU") if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
SET_SOURCE_FILES_PROPERTIES( "lattice.f90" PROPERTIES SET_SOURCE_FILES_PROPERTIES( "lattice.f90" PROPERTIES
COMPILE_FLAGS "-ffree-line-length-240") COMPILE_FLAGS "-ffree-line-length-240")
# long lines for interaction matrix # long lines for interaction matrix
@ -18,10 +18,12 @@ add_library(PREC OBJECT "prec.f90")
add_dependencies(PREC SYSTEM_ROUTINES) add_dependencies(PREC SYSTEM_ROUTINES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>)
if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral") if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_library(DAMASK_INTERFACE OBJECT "spectral_interface.f90") add_library(DAMASK_INTERFACE OBJECT "spectral_interface.f90")
elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM") elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(DAMASK_INTERFACE OBJECT "FEM_interface.f90") add_library(DAMASK_INTERFACE OBJECT "FEM_interface.f90")
else ()
message (FATAL_ERROR "Build target (PROJECT_NAME) is not defined")
endif() endif()
add_dependencies(DAMASK_INTERFACE PREC) add_dependencies(DAMASK_INTERFACE PREC)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>)
@ -47,11 +49,11 @@ add_dependencies(DAMASK_MATH FEsolving)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_MATH>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_MATH>)
# SPECTRAL solver and FEM solver use different mesh files # SPECTRAL solver and FEM solver use different mesh files
if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral") if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_library(MESH OBJECT "mesh.f90") add_library(MESH OBJECT "mesh.f90")
add_dependencies(MESH DAMASK_MATH) add_dependencies(MESH DAMASK_MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM") elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEZoo OBJECT "FEZoo.f90") add_library(FEZoo OBJECT "FEZoo.f90")
add_dependencies(FEZoo DAMASK_MATH) add_dependencies(FEZoo DAMASK_MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
@ -158,7 +160,7 @@ add_library(DAMASK_CPFE OBJECT "CPFEM2.f90")
add_dependencies(DAMASK_CPFE DAMASK_ENGINE) add_dependencies(DAMASK_CPFE DAMASK_ENGINE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_CPFE>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_CPFE>)
if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral") if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_library(SPECTRAL_UTILITIES OBJECT "spectral_utilities.f90") add_library(SPECTRAL_UTILITIES OBJECT "spectral_utilities.f90")
add_dependencies(SPECTRAL_UTILITIES DAMASK_CPFE) add_dependencies(SPECTRAL_UTILITIES DAMASK_CPFE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_UTILITIES>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_UTILITIES>)
@ -170,13 +172,13 @@ if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral")
"spectral_mech_Basic.f90") "spectral_mech_Basic.f90")
add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES) add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>)
if(NOT "${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY") if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES}) add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES})
else() else()
add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90") add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90")
endif() endif()
add_dependencies(DAMASK_spectral SPECTRAL_SOLVER) add_dependencies(DAMASK_spectral SPECTRAL_SOLVER)
elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM") elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEM_UTILITIES OBJECT "FEM_utilities.f90") add_library(FEM_UTILITIES OBJECT "FEM_utilities.f90")
add_dependencies(FEM_UTILITIES DAMASK_CPFE) add_dependencies(FEM_UTILITIES DAMASK_CPFE)

13
src/C_routines.c
View File

@ -11,9 +11,9 @@
int isdirectory_c(const char *dir){ int isdirectory_c(const char *dir){
struct stat statbuf; struct stat statbuf;
if(stat(dir, &statbuf) != 0) if(stat(dir, &statbuf) != 0) /* error */
return 0; return 0; /* return "NO, this is not a directory" */
return S_ISDIR(statbuf.st_mode); return S_ISDIR(statbuf.st_mode); /* 1 => is directory, 0 => this is NOT a directory */
} }
@ -29,7 +29,7 @@ void getcurrentworkdir_c(char cwd[], int *stat ){
} }
void gethostname_c(char hostname[], int *stat ){ void gethostname_c(char hostname[], int *stat){
char hostname_tmp[1024]; char hostname_tmp[1024];
if(gethostname(hostname_tmp, sizeof(hostname_tmp)) == 0){ if(gethostname(hostname_tmp, sizeof(hostname_tmp)) == 0){
strcpy(hostname,hostname_tmp); strcpy(hostname,hostname_tmp);
@ -39,3 +39,8 @@ void gethostname_c(char hostname[], int *stat ){
*stat = 1; *stat = 1;
} }
} }
int chdir_c(const char *dir){
return chdir(dir);
}

2
src/DAMASK_spectral.f90
View File

@ -151,7 +151,7 @@ program DAMASK_spectral
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt) call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>'
write(6,'(/,a)') ' Roters et al., Computational Materials Science, 2018' write(6,'(/,a,/)') ' Roters et al., Computational Materials Science, 2018'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"

133
src/math.f90
View File

@ -1820,6 +1820,8 @@ function math_sampleFiberOri(alpha,beta,FWHM)
integer(pInt):: j,& !< index of smallest component integer(pInt):: j,& !< index of smallest component
i i
allocate(a(0))
allocate(idx(0))
fInC = [sin(alpha(1))*cos(alpha(2)), sin(alpha(1))*sin(alpha(2)), cos(alpha(1))] fInC = [sin(alpha(1))*cos(alpha(2)), sin(alpha(1))*sin(alpha(2)), cos(alpha(1))]
fInS = [sin(beta(1))*cos(beta(2)), sin(beta(1))*sin(beta(2)), cos(beta(1))] fInS = [sin(beta(1))*cos(beta(2)), sin(beta(1))*sin(beta(2)), cos(beta(1))]
@ -2632,135 +2634,4 @@ real(pReal) pure function math_limit(a, left, right)
end function math_limit end function math_limit
!--------------------------------------------------------------------------------------------------
!> @brief Modified Bessel I function of order 0
!> @author John Burkardt
!> @details original version available on https://people.sc.fsu.edu/~jburkardt/f_src/toms715/toms715.html
!--------------------------------------------------------------------------------------------------
real(pReal) function bessel_i0 (x)
use, intrinsic :: IEEE_ARITHMETIC
implicit none
real(pReal), intent(in) :: x
integer(pInt) :: i
real(pReal) :: sump_p, sump_q, xAbs, xx
real(pReal), parameter, dimension(15) :: p_small = real( &
[-5.2487866627945699800e-18, -1.5982226675653184646e-14, -2.6843448573468483278e-11, &
-3.0517226450451067446e-08, -2.5172644670688975051e-05, -1.5453977791786851041e-02, &
-7.0935347449210549190e+00, -2.4125195876041896775e+03, -5.9545626019847898221e+05, &
-1.0313066708737980747e+08, -1.1912746104985237192e+10, -8.4925101247114157499e+11, &
-3.2940087627407749166e+13, -5.5050369673018427753e+14, -2.2335582639474375249e+15], pReal)
real(pReal), parameter, dimension(5) :: q_small = real( &
[-3.7277560179962773046e+03, 6.5158506418655165707e+06, -6.5626560740833869295e+09, &
3.7604188704092954661e+12, -9.7087946179594019126e+14], pReal)
real(pReal), parameter, dimension(8) :: p_large = real( &
[-3.9843750000000000000e-01, 2.9205384596336793945e+00, -2.4708469169133954315e+00, &
4.7914889422856814203e-01, -3.7384991926068969150e-03, -2.6801520353328635310e-03, &
9.9168777670983678974e-05, -2.1877128189032726730e-06], pReal)
real(pReal), parameter, dimension(7) :: q_large = real( &
[-3.1446690275135491500e+01, 8.5539563258012929600e+01, -6.0228002066743340583e+01, &
1.3982595353892851542e+01, -1.1151759188741312645e+00, 3.2547697594819615062e-02, &
-5.5194330231005480228e-04], pReal)
xAbs = abs(x)
argRange: if (xAbs < 5.55e-17_pReal) then
bessel_i0 = 1.0_pReal
else if (xAbs < 15.0_pReal) then argRange
xx = xAbs**2.0_pReal
sump_p = p_small(1)
do i = 2, 15
sump_p = sump_p * xx + p_small(i)
end do
xx = xx - 225.0_pReal
sump_q = ((((xx+q_small(1))*xx+q_small(2))*xx+q_small(3))*xx+q_small(4))*xx+q_small(5)
bessel_i0 = sump_p / sump_q
else if (xAbs <= 713.986_pReal) then argRange
xx = 1.0_pReal / xAbs - 2.0_pReal/30.0_pReal
sump_p = ((((((p_large(1)*xx+p_large(2))*xx+p_large(3))*xx+p_large(4))*xx+ &
p_large(5))*xx+p_large(6))*xx+p_large(7))*xx+p_large(8)
sump_q = ((((((xx+q_large(1))*xx+q_large(2))*xx+q_large(3))*xx+ &
q_large(4))*xx+q_large(5))*xx+q_large(6))*xx+q_large(7)
bessel_i0 = sump_p / sump_q
avoidOverflow: if (xAbs > 698.986_pReal) then
bessel_i0 = ((bessel_i0*exp(xAbs-40.0_pReal)-p_large(1)*exp(xAbs-40.0_pReal))/sqrt(xAbs))*exp(40.0)
else avoidOverflow
bessel_i0 = ((bessel_i0*exp(xAbs)-p_large(1)*exp(xAbs))/sqrt(xAbs))
endif avoidOverflow
else argRange
bessel_i0 = IEEE_value(bessel_i0,IEEE_positive_inf)
end if argRange
end function bessel_i0
!--------------------------------------------------------------------------------------------------
!> @brief Modified Bessel I function of order 1
!> @author John Burkardt
!> @details original version available on https://people.sc.fsu.edu/~jburkardt/f_src/toms715/toms715.html
!--------------------------------------------------------------------------------------------------
real(pReal) function bessel_i1 (x)
use, intrinsic :: IEEE_ARITHMETIC
implicit none
real(pReal), intent(in) :: x
integer(pInt) :: i
real(pReal) :: sump_p, sump_q, xAbs, xx
real(pReal), dimension(15), parameter :: p_small = real( &
[-1.9705291802535139930e-19, -6.5245515583151902910e-16, -1.1928788903603238754e-12, &
-1.4831904935994647675e-09, -1.3466829827635152875e-06, -9.1746443287817501309e-04, &
-4.7207090827310162436e-01, -1.8225946631657315931e+02, -5.1894091982308017540e+04, &
-1.0588550724769347106e+07, -1.4828267606612366099e+09, -1.3357437682275493024e+11, &
-6.9876779648010090070e+12, -1.7732037840791591320e+14, -1.4577180278143463643e+15], pReal)
real(pReal), dimension(5), parameter :: q_small = real( &
[-4.0076864679904189921e+03, 7.4810580356655069138e+06, -8.0059518998619764991e+09, &
4.8544714258273622913e+12, -1.3218168307321442305e+15], pReal)
real(pReal), dimension(8), parameter :: p_large = real( &
[-6.0437159056137600000e-02, 4.5748122901933459000e-01, -4.2843766903304806403e-01, &
9.7356000150886612134e-02, -3.2457723974465568321e-03, -3.6395264712121795296e-04, &
1.6258661867440836395e-05, -3.6347578404608223492e-07], pReal)
real(pReal), dimension(6), parameter :: q_large = real( &
[-3.8806586721556593450e+00, 3.2593714889036996297e+00, -8.5017476463217924408e-01, &
7.4212010813186530069e-02, -2.2835624489492512649e-03, 3.7510433111922824643e-05], pReal)
real(pReal), parameter :: pbar = 3.98437500e-01
xAbs = abs(x)
argRange: if (xAbs < 5.55e-17_pReal) then
bessel_i1 = 0.5_pReal * xAbs
else if (xAbs < 15.0_pReal) then argRange
xx = xAbs**2.0_pReal
sump_p = p_small(1)
do i = 2, 15
sump_p = sump_p * xx + p_small(i)
end do
xx = xx - 225.0_pReal
sump_q = ((((xx+q_small(1))*xx+q_small(2))*xx+q_small(3))*xx+q_small(4)) * xx + q_small(5)
bessel_i1 = (sump_p / sump_q) * xAbs
else if (xAbs <= 713.986_pReal) then argRange
xx = 1.0_pReal / xAbs - 2.0_pReal/30.0_pReal
sump_p = ((((((p_large(1)*xx+p_large(2))*xx+p_large(3))*xx+p_large(4))*xx+&
p_large(5))*xx+p_large(6))*xx+p_large(7))*xx+p_large(8)
sump_q = (((((xx+q_large(1))*xx+q_large(2))*xx+q_large(3))*xx+ q_large(4))*xx+q_large(5))*xx+q_large(6)
bessel_i1 = sump_p / sump_q
avoidOverflow: if (xAbs > 698.986_pReal) then
bessel_i1 = ((bessel_i1 * exp(xAbs-40.0_pReal) + pbar * exp(xAbs-40.0_pReal)) / sqrt(xAbs)) * exp(40.0_pReal)
else avoidOverflow
bessel_i1 = ((bessel_i1 * exp(xAbs) + pbar * exp(xAbs)) / sqrt(xAbs))
endif avoidOverflow
else argRange
bessel_i1 = IEEE_value(bessel_i1,IEEE_positive_inf)
end if argRange
if (x < 0.0_pReal) bessel_i1 = -bessel_i1
end function bessel_i1
end module math end module math

2
src/numerics.f90
View File

@ -108,7 +108,7 @@ module numerics
character(len=64), private :: & character(len=64), private :: &
fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag
character(len=64), protected, public :: & character(len=64), protected, public :: &
spectral_solver = 'basicpetsc' , & !< spectral solution method spectral_solver = 'basic', & !< spectral solution method
spectral_derivative = 'continuous' !< spectral spatial derivative method spectral_derivative = 'continuous' !< spectral spatial derivative method
character(len=1024), protected, public :: & character(len=1024), protected, public :: &
petsc_defaultOptions = '-mech_snes_type ngmres & petsc_defaultOptions = '-mech_snes_type ngmres &

187
src/plastic_isotropic.f90
View File

@ -53,7 +53,7 @@ module plastic_isotropic
dilatation = .false. dilatation = .false.
end type end type
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable, target, private :: param !< containers of constitutive parameters (len Ninstance)
type, private :: tIsotropicState !< internal state aliases type, private :: tIsotropicState !< internal state aliases
real(pReal), pointer, dimension(:) :: & ! scalars along NipcMyInstance real(pReal), pointer, dimension(:) :: & ! scalars along NipcMyInstance
@ -61,20 +61,10 @@ module plastic_isotropic
accumulatedShear accumulatedShear
end type end type
type, private :: tIsotropicAbsTol !< internal alias for abs tolerance in state
real(pReal), pointer :: & ! scalars
flowstress, &
accumulatedShear
end type
type(tIsotropicState), allocatable, dimension(:), private :: & !< state aliases per instance type(tIsotropicState), allocatable, dimension(:), private :: & !< state aliases per instance
state, & state, &
state0, &
dotState dotState
type(tIsotropicAbsTol), allocatable, dimension(:), private :: & !< state aliases per instance
stateAbsTol
public :: & public :: &
plastic_isotropic_init, & plastic_isotropic_init, &
plastic_isotropic_LpAndItsTangent, & plastic_isotropic_LpAndItsTangent, &
@ -130,6 +120,7 @@ subroutine plastic_isotropic_init(fileUnit)
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
type(tParameters), pointer :: p
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: & integer(pInt) :: &
@ -150,8 +141,6 @@ subroutine plastic_isotropic_init(fileUnit)
integer(pInt) :: NipcMyPhase integer(pInt) :: NipcMyPhase
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>' write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'
write(6,'(/,a)') ' Ma et al., Computational Materials Science, 109:323–329, 2015'
write(6,'(/,a)') ' https://doi.org/10.1016/j.commatsci.2015.07.041'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
@ -184,14 +173,11 @@ subroutine plastic_isotropic_init(fileUnit)
endif endif
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
phase = phase + 1_pInt ! advance section counter phase = phase + 1_pInt ! advance section counter
if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then
instance = phase_plasticityInstance(phase) ! count instances of my constitutive law
allocate(param(instance)%outputID(phase_Noutput(phase))) ! allocate space for IDs of every requested output
endif
cycle ! skip to next line cycle ! skip to next line
endif endif
if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
p => param(instance) ! shorthand pointer to parameter object of my constitutive law
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
@ -201,58 +187,58 @@ subroutine plastic_isotropic_init(fileUnit)
select case(outputtag) select case(outputtag)
case ('flowstress') case ('flowstress')
plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
param(instance)%outputID (plastic_isotropic_Noutput(instance)) = flowstress_ID
plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag
p%outputID = [p%outputID,flowstress_ID]
case ('strainrate') case ('strainrate')
plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
param(instance)%outputID (plastic_isotropic_Noutput(instance)) = strainrate_ID
plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag
p%outputID = [p%outputID,strainrate_ID]
end select end select
case ('/dilatation/') case ('/dilatation/')
param(instance)%dilatation = .true. p%dilatation = .true.
case ('tau0') case ('tau0')
param(instance)%tau0 = IO_floatValue(line,chunkPos,2_pInt) p%tau0 = IO_floatValue(line,chunkPos,2_pInt)
case ('gdot0') case ('gdot0')
param(instance)%gdot0 = IO_floatValue(line,chunkPos,2_pInt) p%gdot0 = IO_floatValue(line,chunkPos,2_pInt)
case ('n') case ('n')
param(instance)%n = IO_floatValue(line,chunkPos,2_pInt) p%n = IO_floatValue(line,chunkPos,2_pInt)
case ('h0') case ('h0')
param(instance)%h0 = IO_floatValue(line,chunkPos,2_pInt) p%h0 = IO_floatValue(line,chunkPos,2_pInt)
case ('h0_slope','slopelnrate') case ('h0_slope','slopelnrate')
param(instance)%h0_slopeLnRate = IO_floatValue(line,chunkPos,2_pInt) p%h0_slopeLnRate = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat') case ('tausat')
param(instance)%tausat = IO_floatValue(line,chunkPos,2_pInt) p%tausat = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfita') case ('tausat_sinhfita')
param(instance)%tausat_SinhFitA = IO_floatValue(line,chunkPos,2_pInt) p%tausat_SinhFitA = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitb') case ('tausat_sinhfitb')
param(instance)%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt) p%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitc') case ('tausat_sinhfitc')
param(instance)%tausat_SinhFitC = IO_floatValue(line,chunkPos,2_pInt) p%tausat_SinhFitC = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitd') case ('tausat_sinhfitd')
param(instance)%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt) p%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt)
case ('a', 'w0') case ('a', 'w0')
param(instance)%a = IO_floatValue(line,chunkPos,2_pInt) p%a = IO_floatValue(line,chunkPos,2_pInt)
case ('taylorfactor') case ('taylorfactor')
param(instance)%fTaylor = IO_floatValue(line,chunkPos,2_pInt) p%fTaylor = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_flowstress') case ('atol_flowstress')
param(instance)%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt) p%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_shear') case ('atol_shear')
param(instance)%aTolShear = IO_floatValue(line,chunkPos,2_pInt) p%aTolShear = IO_floatValue(line,chunkPos,2_pInt)
case default case default
@ -261,25 +247,24 @@ subroutine plastic_isotropic_init(fileUnit)
enddo parsingFile enddo parsingFile
allocate(state(maxNinstance)) ! internal state aliases allocate(state(maxNinstance)) ! internal state aliases
allocate(state0(maxNinstance))
allocate(dotState(maxNinstance)) allocate(dotState(maxNinstance))
allocate(stateAbsTol(maxNinstance))
initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop over every plasticity initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop over every plasticity
myPhase: if (phase_plasticity(phase) == PLASTICITY_isotropic_ID) then ! isolate instances of own constitutive description myPhase: if (phase_plasticity(phase) == PLASTICITY_isotropic_ID) then ! isolate instances of own constitutive description
NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc) NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc)
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
p => param(instance)
extmsg = '' extmsg = ''
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity checks
if (param(instance)%aTolShear <= 0.0_pReal) param(instance)%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6 if (p%aTolShear <= 0.0_pReal) p%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6
if (param(instance)%tau0 < 0.0_pReal) extmsg = trim(extmsg)//' tau0' if (p%tau0 < 0.0_pReal) extmsg = trim(extmsg)//' tau0'
if (param(instance)%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' if (p%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
if (param(instance)%n <= 0.0_pReal) extmsg = trim(extmsg)//' n' if (p%n <= 0.0_pReal) extmsg = trim(extmsg)//' n'
if (param(instance)%tausat <= 0.0_pReal) extmsg = trim(extmsg)//' tausat' if (p%tausat <= 0.0_pReal) extmsg = trim(extmsg)//' tausat'
if (param(instance)%a <= 0.0_pReal) extmsg = trim(extmsg)//' a' if (p%a <= 0.0_pReal) extmsg = trim(extmsg)//' a'
if (param(instance)%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//' taylorfactor' if (p%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//' taylorfactor'
if (param(instance)%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_flowstress' if (p%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_flowstress'
if (extmsg /= '') then if (extmsg /= '') then
extmsg = trim(extmsg)//' ('//PLASTICITY_ISOTROPIC_label//')' ! prepare error message identifier extmsg = trim(extmsg)//' ('//PLASTICITY_ISOTROPIC_label//')' ! prepare error message identifier
call IO_error(211_pInt,ip=instance,ext_msg=extmsg) call IO_error(211_pInt,ip=instance,ext_msg=extmsg)
@ -287,7 +272,7 @@ subroutine plastic_isotropic_init(fileUnit)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Determine size of postResults array ! Determine size of postResults array
outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance) outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
select case(param(instance)%outputID(o)) select case(p%outputID(o))
case(flowstress_ID,strainrate_ID) case(flowstress_ID,strainrate_ID)
mySize = 1_pInt mySize = 1_pInt
case default case default
@ -302,7 +287,7 @@ subroutine plastic_isotropic_init(fileUnit)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
sizeDotState = 2_pInt ! flowstress, accumulated_shear sizeDotState = size(["flowstress ","accumulated_shear"])
sizeDeltaState = 0_pInt ! no sudden jumps in state sizeDeltaState = 0_pInt ! no sudden jumps in state
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeState = sizeState
@ -331,31 +316,20 @@ subroutine plastic_isotropic_init(fileUnit)
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! globally required state aliases ! locally defined state aliases and initialization of state0 and aTolState
plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase)
plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases
state(instance)%flowstress => plasticState(phase)%state (1,1:NipcMyPhase) state(instance)%flowstress => plasticState(phase)%state (1,1:NipcMyPhase)
state0(instance)%flowstress => plasticState(phase)%state0 (1,1:NipcMyPhase)
dotState(instance)%flowstress => plasticState(phase)%dotState (1,1:NipcMyPhase) dotState(instance)%flowstress => plasticState(phase)%dotState (1,1:NipcMyPhase)
stateAbsTol(instance)%flowstress => plasticState(phase)%aTolState(1) plasticState(phase)%state0(1,1:NipcMyPhase) = p%tau0
plasticState(phase)%aTolState(1) = p%aTolFlowstress
state(instance)%accumulatedShear => plasticState(phase)%state (2,1:NipcMyPhase) state(instance)%accumulatedShear => plasticState(phase)%state (2,1:NipcMyPhase)
state0(instance)%accumulatedShear => plasticState(phase)%state0 (2,1:NipcMyPhase)
dotState(instance)%accumulatedShear => plasticState(phase)%dotState (2,1:NipcMyPhase) dotState(instance)%accumulatedShear => plasticState(phase)%dotState (2,1:NipcMyPhase)
stateAbsTol(instance)%accumulatedShear => plasticState(phase)%aTolState(2) plasticState(phase)%state0 (2,1:NipcMyPhase) = 0.0_pReal
plasticState(phase)%aTolState(2) = p%aTolShear
!-------------------------------------------------------------------------------------------------- ! global alias
! init state plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase)
state0(instance)%flowstress = param(instance)%tau0 plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase)
state0(instance)%accumulatedShear = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! init absolute state tolerances
stateAbsTol(instance)%flowstress = param(instance)%aTolFlowstress
stateAbsTol(instance)%accumulatedShear = param(instance)%aTolShear
endif myPhase endif myPhase
enddo initializeInstances enddo initializeInstances
@ -400,6 +374,8 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(tParameters), pointer :: p
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
@ -414,7 +390,8 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33) squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) norm_Tstar_dev = sqrt(squarenorm_Tstar_dev)
@ -423,11 +400,11 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dTstar99 = 0.0_pReal dLp_dTstar99 = 0.0_pReal
else else
gamma_dot = param(instance)%gdot0 & gamma_dot = p%gdot0 &
* ( sqrt(1.5_pReal) * norm_Tstar_dev / param(instance)%fTaylor / state(instance)%flowstress(of) ) & * ( sqrt(1.5_pReal) * norm_Tstar_dev / p%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n **p%n
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/param(instance)%fTaylor Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/p%fTaylor
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
@ -441,13 +418,13 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp ! Calculation of the tangent of Lp
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * & dLp_dTstar_3333(k,l,m,n) = (p%n-1.0_pReal) * &
Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
dLp_dTstar99 = math_Plain3333to99(gamma_dot / param(instance)%fTaylor * & dLp_dTstar99 = math_Plain3333to99(gamma_dot / p%fTaylor * &
dLp_dTstar_3333 / norm_Tstar_dev) dLp_dTstar_3333 / norm_Tstar_dev)
end if end if
end subroutine plastic_isotropic_LpAndItsTangent end subroutine plastic_isotropic_LpAndItsTangent
@ -479,9 +456,11 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(tParameters), pointer :: p
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal) :: & real(pReal) :: &
gamma_dot, & !< strainrate gamma_dot, & !< strainrate
norm_Tstar_sph, & !< euclidean norm of Tstar_sph norm_Tstar_sph, & !< euclidean norm of Tstar_sph
squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph
@ -491,34 +470,34 @@ real(pReal) :: &
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33) squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33)
norm_Tstar_sph = sqrt(squarenorm_Tstar_sph) norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
if (param(instance)%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero if (p%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero
gamma_dot = param(instance)%gdot0 & gamma_dot = p%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) & * (sqrt(1.5_pReal) * norm_Tstar_sph / p%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n **p%n
Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/p%fTaylor
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Li ! Calculation of the tangent of Li
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * & dLi_dTstar_3333(k,l,m,n) = (p%n-1.0_pReal) * &
Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * & dLi_dTstar_3333 = gamma_dot / p%fTaylor * &
dLi_dTstar_3333 / norm_Tstar_sph dLi_dTstar_3333 / norm_Tstar_sph
else else
Li = 0.0_pReal Li = 0.0_pReal
dLi_dTstar_3333 = 0.0_pReal dLi_dTstar_3333 = 0.0_pReal
endif endif
end subroutine plastic_isotropic_LiAndItsTangent
end subroutine plastic_isotropic_LiAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -541,6 +520,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(tParameters), pointer :: p
real(pReal), dimension(6) :: & real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: & real(pReal) :: &
@ -554,10 +534,11 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress ! norm of (deviatoric) 2nd Piola-Kirchhoff stress
if (param(instance)%dilatation) then if (p%dilatation) then
norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v)) norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
else else
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
@ -566,26 +547,26 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
end if end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! strain rate ! strain rate
gamma_dot = param(instance)%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & gamma_dot = p%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
/ &!----------------------------------------------------------------------------------- / &!-----------------------------------------------------------------------------------
(param(instance)%fTaylor*state(instance)%flowstress(of) ))**param(instance)%n (p%fTaylor*state(instance)%flowstress(of) ))**p%n
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! hardening coefficient ! hardening coefficient
if (abs(gamma_dot) > 1e-12_pReal) then if (abs(gamma_dot) > 1e-12_pReal) then
if (dEq0(param(instance)%tausat_SinhFitA)) then if (dEq0(p%tausat_SinhFitA)) then
saturation = param(instance)%tausat saturation = p%tausat
else else
saturation = param(instance)%tausat & saturation = p%tausat &
+ asinh( (gamma_dot / param(instance)%tausat_SinhFitA& + asinh( (gamma_dot / p%tausat_SinhFitA&
)**(1.0_pReal / param(instance)%tausat_SinhFitD)& )**(1.0_pReal / p%tausat_SinhFitD)&
)**(1.0_pReal / param(instance)%tausat_SinhFitC) & )**(1.0_pReal / p%tausat_SinhFitC) &
/ ( param(instance)%tausat_SinhFitB & / ( p%tausat_SinhFitB &
* (gamma_dot / param(instance)%gdot0)**(1.0_pReal / param(instance)%n) & * (gamma_dot / p%gdot0)**(1.0_pReal / p%n) &
) )
endif endif
hardening = ( param(instance)%h0 + param(instance)%h0_slopeLnRate * log(gamma_dot) ) & hardening = ( p%h0 + p%h0_slopeLnRate * log(gamma_dot) ) &
* abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**param(instance)%a & * abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**p%a &
* sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation) * sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation)
else else
hardening = 0.0_pReal hardening = 0.0_pReal
@ -614,6 +595,9 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(tParameters), pointer :: p
real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
plastic_isotropic_postResults plastic_isotropic_postResults
@ -629,10 +613,11 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress ! norm of (deviatoric) 2nd Piola-Kirchhoff stress
if (param(instance)%dilatation) then if (p%dilatation) then
norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v)) norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
else else
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
@ -644,15 +629,15 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
plastic_isotropic_postResults = 0.0_pReal plastic_isotropic_postResults = 0.0_pReal
outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance) outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
select case(param(instance)%outputID(o)) select case(p%outputID(o))
case (flowstress_ID) case (flowstress_ID)
plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of) plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of)
c = c + 1_pInt c = c + 1_pInt
case (strainrate_ID) case (strainrate_ID)
plastic_isotropic_postResults(c+1_pInt) = & plastic_isotropic_postResults(c+1_pInt) = &
param(instance)%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & p%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
/ &!---------------------------------------------------------------------------------- / &!----------------------------------------------------------------------------------
(param(instance)%fTaylor * state(instance)%flowstress(of)) ) ** param(instance)%n (p%fTaylor * state(instance)%flowstress(of)) ) ** p%n
c = c + 1_pInt c = c + 1_pInt
end select end select
enddo outputsLoop enddo outputsLoop

106
src/plastic_kinematichardening.f90
View File

@ -12,32 +12,32 @@ module plastic_kinehardening
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_kinehardening_sizePostResults !< cumulative size of post results plastic_kinehardening_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
plastic_kinehardening_sizePostResult !< size of each post result output plastic_kinehardening_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_kinehardening_output !< name of each post result output plastic_kinehardening_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: & integer(pInt), dimension(:), allocatable, target, public :: &
plastic_kinehardening_Noutput !< number of outputs per instance plastic_kinehardening_Noutput !< number of outputs per instance
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_kinehardening_totalNslip !< no. of slip system used in simulation plastic_kinehardening_totalNslip !< no. of slip system used in simulation
integer(pInt), dimension(:,:), allocatable, private :: & integer(pInt), dimension(:,:), allocatable, private :: &
plastic_kinehardening_Nslip !< active number of slip systems per family (input parameter, per family) plastic_kinehardening_Nslip !< active number of slip systems per family (input parameter, per family)
enum, bind(c) enum, bind(c)
enumerator :: undefined_ID, & enumerator :: undefined_ID, &
crss_ID, & !< critical resolved stress crss_ID, & !< critical resolved stress
crss_back_ID, & !< critical resolved back stress crss_back_ID, & !< critical resolved back stress
sense_ID, & !< sense of acting shear stress (-1 or +1) sense_ID, & !< sense of acting shear stress (-1 or +1)
chi0_ID, & !< backstress at last switch of stress sense (positive?) chi0_ID, & !< backstress at last switch of stress sense (positive?)
gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?) gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?)
accshear_ID, & accshear_ID, &
sumGamma_ID, & sumGamma_ID, &
shearrate_ID, & shearrate_ID, &
@ -46,26 +46,26 @@ module plastic_kinehardening
end enum end enum
type, private :: tParameters !< container type for internal constitutive parameters type, private :: tParameters !< container type for internal constitutive parameters
integer(kind(undefined_ID)), dimension(:), allocatable, private :: & integer(kind(undefined_ID)), dimension(:), allocatable, private :: &
outputID !< ID of each post result output outputID !< ID of each post result output
real(pReal) :: & real(pReal) :: &
gdot0, & !< reference shear strain rate for slip (input parameter) gdot0, & !< reference shear strain rate for slip (input parameter)
n_slip, & !< stress exponent for slip (input parameter) n_slip, & !< stress exponent for slip (input parameter)
aTolResistance, & aTolResistance, &
aTolShear aTolShear
real(pReal), dimension(:), allocatable, private :: & real(pReal), dimension(:), allocatable, private :: &
crss0, & !< initial critical shear stress for slip (input parameter, per family) crss0, & !< initial critical shear stress for slip (input parameter, per family)
theta0, & !< initial hardening rate of forward stress for each slip theta0, & !< initial hardening rate of forward stress for each slip
theta1, & !< asymptotic hardening rate of forward stress for each slip > theta1, & !< asymptotic hardening rate of forward stress for each slip >
theta0_b, & !< initial hardening rate of back stress for each slip > theta0_b, & !< initial hardening rate of back stress for each slip >
theta1_b, & !< asymptotic hardening rate of back stress for each slip > theta1_b, & !< asymptotic hardening rate of back stress for each slip >
tau1, & tau1, &
tau1_b, & tau1_b, &
interaction_slipslip, & !< latent hardening matrix interaction_slipslip, & !< latent hardening matrix
nonSchmidCoeff nonSchmidCoeff
real(pReal), dimension(:,:), allocatable, private :: & real(pReal), dimension(:,:), allocatable, private :: &
@ -73,20 +73,20 @@ module plastic_kinehardening
end type end type
type, private :: tKinehardeningState type, private :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress crss, & !< critical resolved stress
crss_back, & !< critical resolved back stress crss_back, & !< critical resolved back stress
sense, & !< sense of acting shear stress (-1 or +1) sense, & !< sense of acting shear stress (-1 or +1)
chi0, & !< backstress at last switch of stress sense chi0, & !< backstress at last switch of stress sense
gamma0, & !< accumulated shear at last switch of stress sense gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear accshear !< accumulated (absolute) shear
real(pReal), pointer, dimension(:) :: & !< scalars along NipcMyInstance real(pReal), pointer, dimension(:) :: & !< scalars along NipcMyInstance
sumGamma !< accumulated shear across all systems sumGamma !< accumulated shear across all systems
end type end type
type(tParameters), dimension(:), allocatable, private :: & type(tParameters), dimension(:), allocatable, private :: &
param !< containers of constitutive parameters (len Ninstance) param !< containers of constitutive parameters (len Ninstance)
type(tKinehardeningState), allocatable, dimension(:), private :: & type(tKinehardeningState), allocatable, dimension(:), private :: &
dotState, & dotState, &
@ -155,9 +155,10 @@ subroutine plastic_kinehardening_init(fileUnit)
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos
integer(kind(undefined_ID)) :: &
output_ID
integer(pInt) :: & integer(pInt) :: &
o, j, k, f, & o, j, k, f, &
output_ID, &
phase, & phase, &
instance, & instance, &
maxNinstance, & maxNinstance, &
@ -177,8 +178,6 @@ subroutine plastic_kinehardening_init(fileUnit)
tag = '', & tag = '', &
line = '', & line = '', &
extmsg = '' extmsg = ''
character(len=64) :: &
outputtag = ''
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_KINEHARDENING_label//' init -+>>>' write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -220,7 +219,6 @@ subroutine plastic_kinehardening_init(fileUnit)
Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) ! maximum number of slip families according to lattice type of current phase Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) ! maximum number of slip families according to lattice type of current phase
Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase)) Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase))
Nchunks_nonSchmid = lattice_NnonSchmid(phase) Nchunks_nonSchmid = lattice_NnonSchmid(phase)
allocate(param(instance)%outputID(phase_Noutput(phase)), source=undefined_ID) ! allocate space for IDs of every requested output
allocate(param(instance)%crss0 (Nchunks_SlipFamilies), source=0.0_pReal) allocate(param(instance)%crss0 (Nchunks_SlipFamilies), source=0.0_pReal)
allocate(param(instance)%tau1 (Nchunks_SlipFamilies), source=0.0_pReal) allocate(param(instance)%tau1 (Nchunks_SlipFamilies), source=0.0_pReal)
allocate(param(instance)%tau1_b (Nchunks_SlipFamilies), source=0.0_pReal) allocate(param(instance)%tau1_b (Nchunks_SlipFamilies), source=0.0_pReal)
@ -236,39 +234,48 @@ subroutine plastic_kinehardening_init(fileUnit)
cycle ! skip to next line cycle ! skip to next line
endif endif
if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_KINEHARDENING_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_KINEHARDENING_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag) select case(tag)
case ('(output)') case ('(output)')
outputtag = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
output_ID = undefined_ID output_ID = undefined_ID
select case(outputtag) select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('resistance') case ('resistance')
output_ID = crss_ID output_ID = crss_ID
case ('backstress') case ('backstress')
output_ID = crss_back_ID output_ID = crss_back_ID
case ('sense') case ('sense')
output_ID = sense_ID output_ID = sense_ID
case ('chi0') case ('chi0')
output_ID = chi0_ID output_ID = chi0_ID
case ('gamma0') case ('gamma0')
output_ID = gamma0_ID output_ID = gamma0_ID
case ('accumulatedshear') case ('accumulatedshear')
output_ID = accshear_ID output_ID = accshear_ID
case ('totalshear') case ('totalshear')
output_ID = sumGamma_ID output_ID = sumGamma_ID
case ('shearrate') case ('shearrate')
output_ID = shearrate_ID output_ID = shearrate_ID
case ('resolvedstress') case ('resolvedstress')
output_ID = resolvedstress_ID output_ID = resolvedstress_ID
end select end select
if (output_ID /= undefined_ID) then if (output_ID /= undefined_ID) then
plastic_kinehardening_Noutput(instance) = plastic_kinehardening_Noutput(instance) + 1_pInt plastic_kinehardening_Noutput(instance) = plastic_kinehardening_Noutput(instance) + 1_pInt
plastic_kinehardening_output(plastic_kinehardening_Noutput(instance),instance) = outputtag plastic_kinehardening_output(plastic_kinehardening_Noutput(instance),instance) = &
param(instance)%outputID (plastic_kinehardening_Noutput(instance)) = output_ID IO_lc(IO_stringValue(line,chunkPos,2_pInt))
param(instance)%outputID = [param(instance)%outputID, output_ID]
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parameters depending on number of slip families ! parameters depending on number of slip families
case ('nslip') case ('nslip')
@ -619,7 +626,6 @@ subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dTstar99, &
math_transpose33 math_transpose33
use lattice, only: & use lattice, only: &
lattice_Sslip, & !< schmid matrix lattice_Sslip, & !< schmid matrix
lattice_Sslip_v, &
lattice_maxNslipFamily, & lattice_maxNslipFamily, &
lattice_NslipSystem, & lattice_NslipSystem, &
lattice_NnonSchmid lattice_NnonSchmid
@ -739,8 +745,6 @@ subroutine plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(phaseAt(ipc,ip,el)))) :: & real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(phaseAt(ipc,ip,el)))) :: &
gdot_pos,gdot_neg, & gdot_pos,gdot_neg, &
tau_pos,tau_neg, & tau_pos,tau_neg, &
@ -799,14 +803,10 @@ end subroutine plastic_kinehardening_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_dotState(Tstar_v,ipc,ip,el) subroutine plastic_kinehardening_dotState(Tstar_v,ipc,ip,el)
use lattice, only: & use lattice, only: &
lattice_Sslip_v, & lattice_maxNslipFamily
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_NnonSchmid
use material, only: & use material, only: &
material_phase, & material_phase, &
phaseAt, phasememberAt, & phaseAt, phasememberAt, &
plasticState, &
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
@ -819,10 +819,8 @@ subroutine plastic_kinehardening_dotState(Tstar_v,ipc,ip,el)
integer(pInt) :: & integer(pInt) :: &
instance,ph, & instance,ph, &
f,i,j,k, & f,i,j, &
index_myFamily,index_otherFamily, &
nSlip, & nSlip, &
offset_accshear, &
of of
real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
@ -873,14 +871,12 @@ end subroutine plastic_kinehardening_dotState
function plastic_kinehardening_postResults(Tstar_v,ipc,ip,el) function plastic_kinehardening_postResults(Tstar_v,ipc,ip,el)
use material, only: & use material, only: &
material_phase, & material_phase, &
plasticState, &
phaseAt, phasememberAt, & phaseAt, phasememberAt, &
phase_plasticityInstance phase_plasticityInstance
use lattice, only: & use lattice, only: &
lattice_Sslip_v, & lattice_Sslip_v, &
lattice_maxNslipFamily, & lattice_maxNslipFamily, &
lattice_NslipSystem, & lattice_NslipSystem
lattice_NnonSchmid
implicit none implicit none
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
@ -896,7 +892,7 @@ function plastic_kinehardening_postResults(Tstar_v,ipc,ip,el)
integer(pInt) :: & integer(pInt) :: &
instance,ph, of, & instance,ph, of, &
nSlip,& nSlip,&
o,f,i,c,j,k, & o,f,i,c,j,&
index_myFamily index_myFamily
real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &

2
src/spectral_mech_Basic.f90
View File

@ -22,7 +22,7 @@ module spectral_mech_basic
private private
character (len=*), parameter, public :: & character (len=*), parameter, public :: &
DAMASK_spectral_SolverBasicPETSC_label = 'basicpetsc' DAMASK_spectral_SolverBasicPETSC_label = 'basic'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types

35
src/system_routines.f90
View File

@ -10,11 +10,12 @@ module system_routines
public :: & public :: &
isDirectory, & isDirectory, &
getCWD, & getCWD, &
getHostName getHostName, &
setCWD
interface interface
function isDirectory_C(path) BIND(C) function isDirectory_C(path) bind(C)
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_INT, & C_INT, &
C_CHAR C_CHAR
@ -38,6 +39,14 @@ interface
integer(C_INT),intent(out) :: stat integer(C_INT),intent(out) :: stat
end subroutine getHostName_C end subroutine getHostName_C
function chdir_C(path) bind(C)
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR
integer(C_INT) :: chdir_C
character(kind=C_CHAR), dimension(1024), intent(in) :: path ! C string is an array
end function chdir_C
end interface end interface
@ -123,5 +132,27 @@ logical function getHostName(str)
end function getHostName end function getHostName
!--------------------------------------------------------------------------------------------------
!> @brief changes the current working directory
!--------------------------------------------------------------------------------------------------
logical function setCWD(path)
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR, &
C_NULL_CHAR
implicit none
character(len=*), intent(in) :: path
character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string is an array
integer :: i
strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength))
do i=1,len(path) ! copy array components
strFixedLength(i)=path(i:i)
enddo
setCWD=merge(.True.,.False.,chdir_C(strFixedLength) /= 0_C_INT)
end function setCWD
end module system_routines end module system_routines