Merge branch 'development' into 23_BasticPETSc_2_PETSc

This commit is contained in:
Jaeyong Jung 2018-06-04 17:03:13 +02:00
commit 81dffef21c
12 changed files with 230 additions and 336 deletions

View File

@ -5,8 +5,8 @@ cmake_minimum_required (VERSION 2.8.8 FATAL_ERROR)
#---------------------------------------------------------------------------------------
# Find PETSc from system environment
set(PETSC_DIR $ENV{PETSC_DIR})
if ("${PETSC_DIR}" STREQUAL "")
message (FATAL_ERROR "PETSC_DIR is not defined")
if (PETSC_DIR STREQUAL "")
message (FATAL_ERROR "PETSc location (PETSC_DIR) is not defined")
endif ()
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
# DAMASK solver defines project to build
if ("${DAMASK_SOLVER}" STREQUAL "SPECTRAL")
if (DAMASK_SOLVER STREQUAL "SPECTRAL")
project (DAMASK_spectral Fortran C)
add_definitions (-DSpectral)
message ("Building Spectral Solver\n")
elseif ("${DAMASK_SOLVER}" STREQUAL "FEM")
elseif (DAMASK_SOLVER STREQUAL "FEM")
project (DAMASK_FEM Fortran C)
add_definitions (-DFEM)
message ("Building FEM Solver\n")
else ()
message (FATAL_ERROR "Build target (DAMASK_SOLVER) is not defined")
endif ()
# set linker commands (needs to be done after defining the project)
set (CMAKE_LINKER "${PETSC_LINKER}")
if ("${CMAKE_BUILD_TYPE}" STREQUAL "")
if (CMAKE_BUILD_TYPE STREQUAL "")
set (CMAKE_BUILD_TYPE "RELEASE")
endif ()
# 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 (PARALLEL "OFF")
set (OPTI "OFF")
elseif ("${CMAKE_BUILD_TYPE}" STREQUAL "RELEASE")
elseif (CMAKE_BUILD_TYPE STREQUAL "RELEASE")
set (PARALLEL "ON")
set (OPTI "DEFENSIVE")
elseif ("${CMAKE_BUILD_TYPE}" STREQUAL "PERFORMANCE")
elseif (CMAKE_BUILD_TYPE STREQUAL "PERFORMANCE")
set (PARALLEL "ON")
set (OPTI "AGGRESSIVE")
endif ()
# $OPTIMIZATION takes precedence over $BUILD_TYPE defaults
if ("${OPTIMIZATION}" STREQUAL "")
if (OPTIMIZATION STREQUAL "")
set (OPTIMIZATION "${OPTI}")
else ()
set (OPTIMIZATION "${OPTIMIZATION}")
endif ()
# $OPENMP takes precedence over $BUILD_TYPE defaults
if ("${OPENMP}" STREQUAL "")
if (OPENMP STREQUAL "")
set (OPENMP "${PARALLEL}")
else ()
set(OPENMP "${OPENMP}")
endif ()
# 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")
endif ()
@ -188,19 +190,19 @@ set (DAMASK_INCLUDE_FLAGS "${DAMASK_INCLUDE_FLAGS} ${PETSC_INCLUDES}")
###################################################################################################
# Intel Compiler
###################################################################################################
if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel")
if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
if (OPENMP)
set (OPENMP_FLAGS "-qopenmp -parallel")
endif ()
if ("${OPTIMIZATION}" STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0 -no-ip")
elseif ("${OPTIMIZATION}" STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
elseif ("${OPTIMIZATION}" STREQUAL "AGGRESSIVE")
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"
if (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0 -no-ip")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
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"
endif ()
# -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
###################################################################################################
elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU")
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
if (OPENMP)
set (OPENMP_FLAGS "-fopenmp")
endif ()
if ("${OPTIMIZATION}" STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0" )
elseif ("${OPTIMIZATION}" STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
elseif ("${OPTIMIZATION}" STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-O3 -ffast-math -funroll-loops -ftree-vectorize")
if (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0" )
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-O3 -ffast-math -funroll-loops -ftree-vectorize")
endif ()
set (STANDARD_CHECK "-std=f2008ts -pedantic-errors" )
@ -443,12 +445,15 @@ elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU")
# 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)
else ()
message (FATAL_ERROR "Compiler type (CMAKE_Fortran_COMPILER_ID) not recognized")
endif ()
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}")
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_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}")
endif ()
@ -464,15 +469,15 @@ message ("Fortran Linker Command:\n${CMAKE_Fortran_LINK_EXECUTABLE}\n")
add_subdirectory (src)
# INSTALL BUILT BINARIES
if ("${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY")
if (CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
exec_program (mktemp ARGS -d OUTPUT_VARIABLE BLACK_HOLE)
install (PROGRAMS ${PROJECT_BINARY_DIR}/src/prec.mod
DESTINATION ${BLACK_HOLE})
else ()
if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral")
if (PROJECT_NAME STREQUAL "DAMASK_spectral")
install (PROGRAMS ${PROJECT_BINARY_DIR}/src/DAMASK_spectral
DESTINATION ${CMAKE_INSTALL_PREFIX})
elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM")
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
install (PROGRAMS ${PROJECT_BINARY_DIR}/src/DAMASK_FEM
DESTINATION ${CMAKE_INSTALL_PREFIX})
endif ()

View File

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

@ -1 +1 @@
Subproject commit e8eea0593728e1c1646344975dd081fea04842bd
Subproject commit 701d63b0e11a653797afe260d1dfc12e2a390d6f

View File

@ -1 +1 @@
v2.0.1-1227-g7f0275e
v2.0.2-35-g4b5401e

View File

@ -1,5 +1,5 @@
# 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
COMPILE_FLAGS "-ffree-line-length-240")
# long lines for interaction matrix
@ -18,10 +18,12 @@ add_library(PREC OBJECT "prec.f90")
add_dependencies(PREC SYSTEM_ROUTINES)
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")
elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM")
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(DAMASK_INTERFACE OBJECT "FEM_interface.f90")
else ()
message (FATAL_ERROR "Build target (PROJECT_NAME) is not defined")
endif()
add_dependencies(DAMASK_INTERFACE PREC)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>)
@ -47,11 +49,11 @@ add_dependencies(DAMASK_MATH FEsolving)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_MATH>)
# 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_dependencies(MESH DAMASK_MATH)
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_dependencies(FEZoo DAMASK_MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
@ -158,7 +160,7 @@ add_library(DAMASK_CPFE OBJECT "CPFEM2.f90")
add_dependencies(DAMASK_CPFE DAMASK_ENGINE)
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_dependencies(SPECTRAL_UTILITIES DAMASK_CPFE)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_UTILITIES>)
@ -170,13 +172,13 @@ if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral")
"spectral_mech_Basic.f90")
add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES)
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})
else()
add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90")
endif()
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_dependencies(FEM_UTILITIES DAMASK_CPFE)

View File

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

View File

@ -151,7 +151,7 @@ program DAMASK_spectral
! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
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()
#include "compilation_info.f90"

View File

@ -1820,6 +1820,8 @@ function math_sampleFiberOri(alpha,beta,FWHM)
integer(pInt):: j,& !< index of smallest component
i
allocate(a(0))
allocate(idx(0))
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))]
@ -2632,135 +2634,4 @@ real(pReal) pure function math_limit(a, left, right)
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

View File

@ -108,7 +108,7 @@ module numerics
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
character(len=64), protected, public :: &
spectral_solver = 'basic' , & !< spectral solution method
spectral_solver = 'basic', & !< spectral solution method
spectral_derivative = 'continuous' !< spectral spatial derivative method
character(len=1024), protected, public :: &
petsc_defaultOptions = '-mech_snes_type ngmres &

View File

@ -53,7 +53,7 @@ module plastic_isotropic
dilatation = .false.
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
real(pReal), pointer, dimension(:) :: & ! scalars along NipcMyInstance
@ -61,20 +61,10 @@ module plastic_isotropic
accumulatedShear
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
state, &
state0, &
dotState
type(tIsotropicAbsTol), allocatable, dimension(:), private :: & !< state aliases per instance
stateAbsTol
public :: &
plastic_isotropic_init, &
plastic_isotropic_LpAndItsTangent, &
@ -130,6 +120,7 @@ subroutine plastic_isotropic_init(fileUnit)
implicit none
integer(pInt), intent(in) :: fileUnit
type(tParameters), pointer :: p
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: &
@ -150,8 +141,6 @@ subroutine plastic_isotropic_init(fileUnit)
integer(pInt) :: NipcMyPhase
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'
write(6,'(/,a)') ' Ma et al., Computational Materials Science, 109:323329, 2015'
write(6,'(/,a)') ' https://doi.org/10.1016/j.commatsci.2015.07.041'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
@ -184,14 +173,11 @@ subroutine plastic_isotropic_init(fileUnit)
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
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
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
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)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
@ -201,58 +187,58 @@ subroutine plastic_isotropic_init(fileUnit)
select case(outputtag)
case ('flowstress')
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
p%outputID = [p%outputID,flowstress_ID]
case ('strainrate')
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
p%outputID = [p%outputID,strainrate_ID]
end select
case ('/dilatation/')
param(instance)%dilatation = .true.
p%dilatation = .true.
case ('tau0')
param(instance)%tau0 = IO_floatValue(line,chunkPos,2_pInt)
p%tau0 = IO_floatValue(line,chunkPos,2_pInt)
case ('gdot0')
param(instance)%gdot0 = IO_floatValue(line,chunkPos,2_pInt)
p%gdot0 = IO_floatValue(line,chunkPos,2_pInt)
case ('n')
param(instance)%n = IO_floatValue(line,chunkPos,2_pInt)
p%n = IO_floatValue(line,chunkPos,2_pInt)
case ('h0')
param(instance)%h0 = IO_floatValue(line,chunkPos,2_pInt)
p%h0 = IO_floatValue(line,chunkPos,2_pInt)
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')
param(instance)%tausat = IO_floatValue(line,chunkPos,2_pInt)
p%tausat = IO_floatValue(line,chunkPos,2_pInt)
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')
param(instance)%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt)
p%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt)
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')
param(instance)%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt)
p%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt)
case ('a', 'w0')
param(instance)%a = IO_floatValue(line,chunkPos,2_pInt)
p%a = IO_floatValue(line,chunkPos,2_pInt)
case ('taylorfactor')
param(instance)%fTaylor = IO_floatValue(line,chunkPos,2_pInt)
p%fTaylor = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_flowstress')
param(instance)%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt)
p%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_shear')
param(instance)%aTolShear = IO_floatValue(line,chunkPos,2_pInt)
p%aTolShear = IO_floatValue(line,chunkPos,2_pInt)
case default
@ -261,25 +247,24 @@ subroutine plastic_isotropic_init(fileUnit)
enddo parsingFile
allocate(state(maxNinstance)) ! internal state aliases
allocate(state0(maxNinstance))
allocate(dotState(maxNinstance))
allocate(stateAbsTol(maxNinstance))
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
NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc)
instance = phase_plasticityInstance(phase)
p => param(instance)
extmsg = ''
!--------------------------------------------------------------------------------------------------
! sanity checks
if (param(instance)%aTolShear <= 0.0_pReal) param(instance)%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6
if (param(instance)%tau0 < 0.0_pReal) extmsg = trim(extmsg)//' tau0'
if (param(instance)%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
if (param(instance)%n <= 0.0_pReal) extmsg = trim(extmsg)//' n'
if (param(instance)%tausat <= 0.0_pReal) extmsg = trim(extmsg)//' tausat'
if (param(instance)%a <= 0.0_pReal) extmsg = trim(extmsg)//' a'
if (param(instance)%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//' taylorfactor'
if (param(instance)%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_flowstress'
if (p%aTolShear <= 0.0_pReal) p%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6
if (p%tau0 < 0.0_pReal) extmsg = trim(extmsg)//' tau0'
if (p%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
if (p%n <= 0.0_pReal) extmsg = trim(extmsg)//' n'
if (p%tausat <= 0.0_pReal) extmsg = trim(extmsg)//' tausat'
if (p%a <= 0.0_pReal) extmsg = trim(extmsg)//' a'
if (p%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//' taylorfactor'
if (p%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_flowstress'
if (extmsg /= '') then
extmsg = trim(extmsg)//' ('//PLASTICITY_ISOTROPIC_label//')' ! prepare error message identifier
call IO_error(211_pInt,ip=instance,ext_msg=extmsg)
@ -287,7 +272,7 @@ subroutine plastic_isotropic_init(fileUnit)
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
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)
mySize = 1_pInt
case default
@ -302,7 +287,7 @@ subroutine plastic_isotropic_init(fileUnit)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
sizeDotState = 2_pInt ! flowstress, accumulated_shear
sizeDotState = size(["flowstress ","accumulated_shear"])
sizeDeltaState = 0_pInt ! no sudden jumps in state
sizeState = sizeDotState + sizeDeltaState
plasticState(phase)%sizeState = sizeState
@ -331,31 +316,20 @@ subroutine plastic_isotropic_init(fileUnit)
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal)
!--------------------------------------------------------------------------------------------------
! globally required state aliases
plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase)
plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase)
! locally defined state aliases and initialization of state0 and aTolState
!--------------------------------------------------------------------------------------------------
! locally defined state aliases
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)
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)
state0(instance)%accumulatedShear => plasticState(phase)%state0 (2,1:NipcMyPhase)
dotState(instance)%accumulatedShear => plasticState(phase)%dotState (2,1:NipcMyPhase)
stateAbsTol(instance)%accumulatedShear => plasticState(phase)%aTolState(2)
!--------------------------------------------------------------------------------------------------
! init state
state0(instance)%flowstress = param(instance)%tau0
state0(instance)%accumulatedShear = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! init absolute state tolerances
stateAbsTol(instance)%flowstress = param(instance)%aTolFlowstress
stateAbsTol(instance)%accumulatedShear = param(instance)%aTolShear
plasticState(phase)%state0 (2,1:NipcMyPhase) = 0.0_pReal
plasticState(phase)%aTolState(2) = p%aTolShear
! global alias
plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase)
plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase)
endif myPhase
enddo initializeInstances
@ -400,6 +374,8 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
ip, & !< integration point
el !< element
type(tParameters), pointer :: p
real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
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
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
squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
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
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = param(instance)%gdot0 &
* ( sqrt(1.5_pReal) * norm_Tstar_dev / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n
gamma_dot = p%gdot0 &
* ( sqrt(1.5_pReal) * norm_Tstar_dev / p%fTaylor / state(instance)%flowstress(of) ) &
**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 &
.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
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
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
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_dTstar99 = math_Plain3333to99(gamma_dot / param(instance)%fTaylor * &
dLp_dTstar99 = math_Plain3333to99(gamma_dot / p%fTaylor * &
dLp_dTstar_3333 / norm_Tstar_dev)
end if
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
el !< element
type(tParameters), pointer :: p
real(pReal), dimension(3,3) :: &
Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal) :: &
real(pReal) :: &
gamma_dot, & !< strainrate
norm_Tstar_sph, & !< 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
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
squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33)
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
gamma_dot = param(instance)%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n
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 = p%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / p%fTaylor / state(instance)%flowstress(of) ) &
**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
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
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 = gamma_dot / param(instance)%fTaylor * &
dLi_dTstar_3333 = gamma_dot / p%fTaylor * &
dLi_dTstar_3333 / norm_Tstar_sph
else
Li = 0.0_pReal
dLi_dTstar_3333 = 0.0_pReal
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
ip, & !< integration point
el !< element
type(tParameters), pointer :: p
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
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
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
!--------------------------------------------------------------------------------------------------
! 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))
else
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
!--------------------------------------------------------------------------------------------------
! 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
if (abs(gamma_dot) > 1e-12_pReal) then
if (dEq0(param(instance)%tausat_SinhFitA)) then
saturation = param(instance)%tausat
if (dEq0(p%tausat_SinhFitA)) then
saturation = p%tausat
else
saturation = param(instance)%tausat &
+ asinh( (gamma_dot / param(instance)%tausat_SinhFitA&
)**(1.0_pReal / param(instance)%tausat_SinhFitD)&
)**(1.0_pReal / param(instance)%tausat_SinhFitC) &
/ ( param(instance)%tausat_SinhFitB &
* (gamma_dot / param(instance)%gdot0)**(1.0_pReal / param(instance)%n) &
saturation = p%tausat &
+ asinh( (gamma_dot / p%tausat_SinhFitA&
)**(1.0_pReal / p%tausat_SinhFitD)&
)**(1.0_pReal / p%tausat_SinhFitC) &
/ ( p%tausat_SinhFitB &
* (gamma_dot / p%gdot0)**(1.0_pReal / p%n) &
)
endif
hardening = ( param(instance)%h0 + param(instance)%h0_slopeLnRate * log(gamma_dot) ) &
* abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**param(instance)%a &
hardening = ( p%h0 + p%h0_slopeLnRate * log(gamma_dot) ) &
* abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**p%a &
* sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation)
else
hardening = 0.0_pReal
@ -614,6 +595,9 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(tParameters), pointer :: p
real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
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
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
!--------------------------------------------------------------------------------------------------
! 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))
else
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
outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
select case(param(instance)%outputID(o))
select case(p%outputID(o))
case (flowstress_ID)
plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of)
c = c + 1_pInt
case (strainrate_ID)
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
end select
enddo outputsLoop

View File

@ -12,32 +12,32 @@ module plastic_kinehardening
implicit none
private
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 :: &
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 :: &
plastic_kinehardening_output !< name of each post result output
plastic_kinehardening_output !< name of each post result output
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 :: &
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 :: &
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)
enumerator :: undefined_ID, &
crss_ID, & !< critical resolved stress
crss_back_ID, & !< critical resolved back stress
sense_ID, & !< sense of acting shear stress (-1 or +1)
chi0_ID, & !< backstress at last switch of stress sense (positive?)
gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?)
crss_ID, & !< critical resolved stress
crss_back_ID, & !< critical resolved back stress
sense_ID, & !< sense of acting shear stress (-1 or +1)
chi0_ID, & !< backstress at last switch of stress sense (positive?)
gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?)
accshear_ID, &
sumGamma_ID, &
shearrate_ID, &
@ -46,26 +46,26 @@ module plastic_kinehardening
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 :: &
outputID !< ID of each post result output
outputID !< ID of each post result output
real(pReal) :: &
gdot0, & !< reference shear strain rate for slip (input parameter)
n_slip, & !< stress exponent for slip (input parameter)
gdot0, & !< reference shear strain rate for slip (input parameter)
n_slip, & !< stress exponent for slip (input parameter)
aTolResistance, &
aTolShear
real(pReal), dimension(:), allocatable, private :: &
crss0, & !< initial critical shear stress for slip (input parameter, per family)
theta0, & !< initial 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 >
theta1_b, & !< asymptotic hardening rate of back stress for each slip >
crss0, & !< initial critical shear stress for slip (input parameter, per family)
theta0, & !< initial 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 >
theta1_b, & !< asymptotic hardening rate of back stress for each slip >
tau1, &
tau1_b, &
interaction_slipslip, & !< latent hardening matrix
interaction_slipslip, & !< latent hardening matrix
nonSchmidCoeff
real(pReal), dimension(:,:), allocatable, private :: &
@ -73,20 +73,20 @@ module plastic_kinehardening
end type
type, private :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress
crss_back, & !< critical resolved back stress
sense, & !< sense of acting shear stress (-1 or +1)
chi0, & !< backstress at last switch of stress sense
gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
crss, & !< critical resolved stress
crss_back, & !< critical resolved back stress
sense, & !< sense of acting shear stress (-1 or +1)
chi0, & !< backstress at last switch of stress sense
gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear
real(pReal), pointer, dimension(:) :: & !< scalars along NipcMyInstance
sumGamma !< accumulated shear across all systems
real(pReal), pointer, dimension(:) :: & !< scalars along NipcMyInstance
sumGamma !< accumulated shear across all systems
end type
type(tParameters), dimension(:), allocatable, private :: &
param !< containers of constitutive parameters (len Ninstance)
param !< containers of constitutive parameters (len Ninstance)
type(tKinehardeningState), allocatable, dimension(:), private :: &
dotState, &
@ -155,9 +155,10 @@ subroutine plastic_kinehardening_init(fileUnit)
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(kind(undefined_ID)) :: &
output_ID
integer(pInt) :: &
o, j, k, f, &
output_ID, &
phase, &
instance, &
maxNinstance, &
@ -177,8 +178,6 @@ subroutine plastic_kinehardening_init(fileUnit)
tag = '', &
line = '', &
extmsg = ''
character(len=64) :: &
outputtag = ''
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'
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_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,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)%tau1 (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
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
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
outputtag = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
output_ID = undefined_ID
select case(outputtag)
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('resistance')
output_ID = crss_ID
case ('backstress')
output_ID = crss_back_ID
case ('sense')
output_ID = sense_ID
case ('chi0')
output_ID = chi0_ID
case ('gamma0')
output_ID = gamma0_ID
case ('accumulatedshear')
output_ID = accshear_ID
case ('totalshear')
output_ID = sumGamma_ID
case ('shearrate')
output_ID = shearrate_ID
case ('resolvedstress')
output_ID = resolvedstress_ID
end select
if (output_ID /= undefined_ID) then
plastic_kinehardening_Noutput(instance) = plastic_kinehardening_Noutput(instance) + 1_pInt
plastic_kinehardening_output(plastic_kinehardening_Noutput(instance),instance) = outputtag
param(instance)%outputID (plastic_kinehardening_Noutput(instance)) = output_ID
plastic_kinehardening_output(plastic_kinehardening_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
param(instance)%outputID = [param(instance)%outputID, output_ID]
endif
!--------------------------------------------------------------------------------------------------
! parameters depending on number of slip families
case ('nslip')
@ -619,7 +626,6 @@ subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dTstar99, &
math_transpose33
use lattice, only: &
lattice_Sslip, & !< schmid matrix
lattice_Sslip_v, &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_NnonSchmid
@ -739,8 +745,6 @@ subroutine plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el)
ipc, & !< component-ID of integration point
ip, & !< integration point
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)))) :: &
gdot_pos,gdot_neg, &
tau_pos,tau_neg, &
@ -799,14 +803,10 @@ end subroutine plastic_kinehardening_deltaState
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_dotState(Tstar_v,ipc,ip,el)
use lattice, only: &
lattice_Sslip_v, &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_NnonSchmid
lattice_maxNslipFamily
use material, only: &
material_phase, &
phaseAt, phasememberAt, &
plasticState, &
phase_plasticityInstance
implicit none
@ -819,10 +819,8 @@ subroutine plastic_kinehardening_dotState(Tstar_v,ipc,ip,el)
integer(pInt) :: &
instance,ph, &
f,i,j,k, &
index_myFamily,index_otherFamily, &
f,i,j, &
nSlip, &
offset_accshear, &
of
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)
use material, only: &
material_phase, &
plasticState, &
phaseAt, phasememberAt, &
phase_plasticityInstance
use lattice, only: &
lattice_Sslip_v, &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_NnonSchmid
lattice_NslipSystem
implicit none
real(pReal), dimension(6), intent(in) :: &
@ -896,7 +892,7 @@ function plastic_kinehardening_postResults(Tstar_v,ipc,ip,el)
integer(pInt) :: &
instance,ph, of, &
nSlip,&
o,f,i,c,j,k, &
o,f,i,c,j,&
index_myFamily
real(pReal), dimension(plastic_kinehardening_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &

View File

@ -10,11 +10,12 @@ module system_routines
public :: &
isDirectory, &
getCWD, &
getHostName
getHostName, &
setCWD
interface
function isDirectory_C(path) BIND(C)
function isDirectory_C(path) bind(C)
use, intrinsic :: ISO_C_Binding, only: &
C_INT, &
C_CHAR
@ -38,6 +39,14 @@ interface
integer(C_INT),intent(out) :: stat
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
@ -123,5 +132,27 @@ logical function getHostName(str)
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