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

Merge branch 'release_bugfix' into 'development' 

Improve error reports in CMake script

See merge request damask/DAMASK!26
This commit is contained in:
Martin Diehl 2018-06-01 23:34:28 +02:00
parent 4bc89ad5aa 1ad700f7a0
commit 4b5401e10c
5 changed files with 96 additions and 225 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

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)

131
src/math.f90
View File

@ -2632,135 +2632,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

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)))) :: &