diff --git a/CMakeLists.txt b/CMakeLists.txt index cd33a4b31..ee376ef02 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 () diff --git a/DAMASK_env.csh b/DAMASK_env.csh deleted file mode 120000 index e8a0a2c05..000000000 --- a/DAMASK_env.csh +++ /dev/null @@ -1 +0,0 @@ -env/DAMASK.csh \ No newline at end of file diff --git a/DAMASK_env.sh b/DAMASK_env.sh deleted file mode 120000 index 264b07d52..000000000 --- a/DAMASK_env.sh +++ /dev/null @@ -1 +0,0 @@ -env/DAMASK.sh \ No newline at end of file diff --git a/DAMASK_env.zsh b/DAMASK_env.zsh deleted file mode 120000 index cf3a247ef..000000000 --- a/DAMASK_env.zsh +++ /dev/null @@ -1 +0,0 @@ -env/DAMASK.zsh \ No newline at end of file diff --git a/Makefile b/Makefile index 0ba67963e..4149078d1 100644 --- a/Makefile +++ b/Makefile @@ -25,7 +25,6 @@ build/FEM: .PHONY: marc marc: - @./installation/symLink_Code.sh @./installation/mods_MarcMentat/apply_DAMASK_modifications.sh ${MAKEFLAGS} .PHONY: clean diff --git a/PRIVATE b/PRIVATE index 74bae0779..cd02f6c1a 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 74bae0779137ad7ecbab44624b602e32f48a899d +Subproject commit cd02f6c1a481491eb4517651516b8311348b4777 diff --git a/VERSION b/VERSION index f3b15f3f8..dd4a5b4a2 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2 +v2.0.2-37-gc903880 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9789ec67d..9a828669b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 $) -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 $) @@ -47,11 +49,11 @@ add_dependencies(DAMASK_MATH FEsolving) list(APPEND OBJECTFILES $) # 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 $) -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 $) @@ -158,7 +160,7 @@ add_library(DAMASK_CPFE OBJECT "CPFEM2.f90") add_dependencies(DAMASK_CPFE DAMASK_ENGINE) list(APPEND OBJECTFILES $) -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 $) @@ -170,13 +172,13 @@ if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral") "spectral_mech_Basic.f90") add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES) list(APPEND OBJECTFILES $) - 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) diff --git a/src/C_routines.c b/src/C_routines.c index 5bc09745f..e3891765a 100644 --- a/src/C_routines.c +++ b/src/C_routines.c @@ -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); +} diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index ba7c93144..f340b683d 100644 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -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" diff --git a/src/math.f90 b/src/math.f90 index 8f228582b..4b8a00646 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -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 diff --git a/src/numerics.f90 b/src/numerics.f90 index 3a91a36bd..27b04cd67 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -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 = 'basicpetsc' , & !< 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 & diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index a3a1d5caf..ad62ed398 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -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:323–329, 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 diff --git a/src/plastic_kinematichardening.f90 b/src/plastic_kinematichardening.f90 index 434495043..7aba2730b 100644 --- a/src/plastic_kinematichardening.f90 +++ b/src/plastic_kinematichardening.f90 @@ -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)))) :: & diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index 8fd2b7d96..1c949bb7b 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -22,7 +22,7 @@ module spectral_mech_basic private character (len=*), parameter, public :: & - DAMASK_spectral_SolverBasicPETSC_label = 'basicpetsc' + DAMASK_spectral_SolverBasicPETSC_label = 'basic' !-------------------------------------------------------------------------------------------------- ! derived types diff --git a/src/system_routines.f90 b/src/system_routines.f90 index 07e12a20b..2740011b4 100644 --- a/src/system_routines.f90 +++ b/src/system_routines.f90 @@ -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