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

Merge branch '30_parsePhasePartOnce' into 20-NewStyleDislotwin

This commit is contained in:
Sharan Roongta 2018-06-25 15:36:23 +02:00
parent 277cc59069 300b2827b2
commit 57c6f44dfc
58 changed files with 1588 additions and 1400 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
.gitlab-ci.yml
View File

@ -52,23 +52,23 @@ variables:
IntelCompiler16_0: "Compiler/Intel/16.0 Libraries/IMKL/2016"
IntelCompiler17_0: "Compiler/Intel/17.0 Libraries/IMKL/2017"
IntelCompiler18_1: "Compiler/Intel/18.1 Libraries/IMKL/2018"
GNUCompiler5_3: "Compiler/GNU/5.3"
GNUCompiler7_3: "Compiler/GNU/7.3"
# ------------ Defaults ----------------------------------------------
IntelCompiler: "$IntelCompiler18_1"
GNUCompiler: "$GNUCompiler5_3"
GNUCompiler: "$GNUCompiler7_3"
# ++++++++++++ MPI +++++++++++++++++++++++++++++++++++++++++++++++++++
MPICH3_2Intel17_0: "MPI/Intel/17.0/MPICH/3.2"
MPICH3_2Intel18_1: "MPI/Intel/18.1/MPICH/3.2.1"
MPICH3_2GNU5_3: "MPI/GNU/5.3/MPICH/3.2.1"
MPICH3_2GNU7_3: "MPI/GNU/7.3/MPICH/3.2.1"
# ------------ Defaults ----------------------------------------------
MPICH_GNU: "$MPICH3_2GNU5_3"
MPICH_GNU: "$MPICH3_2GNU7_3"
MPICH_Intel: "$MPICH3_2Intel18_1"
# ++++++++++++ PETSc +++++++++++++++++++++++++++++++++++++++++++++++++
PETSc3_9_1MPICH3_2Intel18_1: "Libraries/PETSc/3.9.1/Intel-18.1-MPICH-3.2.1"
PETSc3_9_1MPICH3_2GNU5_3: "Libraries/PETSc/3.9.1/GNU-5.3-MPICH-3.2.1"
PETSc3_9_1MPICH3_2GNU7_3: "Libraries/PETSc/3.9.1/GNU-7.3-MPICH-3.2.1"
# ------------ Defaults ----------------------------------------------
PETSc_MPICH_Intel: "$PETSc3_9_1MPICH3_2Intel18_1"
PETSc_MPICH_GNU: "$PETSc3_9_1MPICH3_2GNU5_3"
PETSc_MPICH_GNU: "$PETSc3_9_1MPICH3_2GNU7_3"
# ++++++++++++ FEM +++++++++++++++++++++++++++++++++++++++++++++++++++
Abaqus2016: "FEM/Abaqus/2016"
Abaqus2017: "FEM/Abaqus/2017"

51
CMakeLists.txt
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,17 +190,17 @@ 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")
if (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0 -no-ip")
elseif ("${OPTIMIZATION}" STREQUAL "DEFENSIVE")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
elseif ("${OPTIMIZATION}" STREQUAL "AGGRESSIVE")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost")
# -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost"
endif ()
@ -308,17 +310,17 @@ 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")
if (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0" )
elseif ("${OPTIMIZATION}" STREQUAL "DEFENSIVE")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
elseif ("${OPTIMIZATION}" STREQUAL "AGGRESSIVE")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-O3 -ffast-math -funroll-loops -ftree-vectorize")
endif ()
@ -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 ()

1
Makefile
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

2
PRIVATE

@ -1 +1 @@
Subproject commit cd02f6c1a481491eb4517651516b8311348b4777
Subproject commit aead92902b3a0cf3404be9c552bfec918d7aaffb

2
VERSION
View File

@ -1 +1 @@
v2.0.2-22-g60e30e4
v2.0.2-48-gaebb06e

2
examples/SpectralMethod/Polycrystal/material.config
View File

@ -3,7 +3,7 @@
#-------------------#
[SX]
type none
mech none
#-------------------#
<crystallite>

56
examples/SpectralMethod/Polycrystal/orientation.seeds
View File

@ -1,56 +0,0 @@
5 header
seeds_fromRandom v2.0.1-1138-gfcac08c -N 50 -g 128 128 128
grid a 128 b 128 c 128
microstructures 50
randomSeed 3336946323
1_pos 2_pos 3_pos 1_euler 2_euler 3_euler microstructure
0.54457843603947365 0.84911587396210719 0.34846714169395199 146.18027121829002 137.38970467457548 64.889274068548971 1.0
0.30082506347847232 0.98313838966599176 0.44557226838658942 277.4997516434205 39.360506400353323 71.246613676352894 2.0
0.40772634005027159 0.9616152434202665 0.058204060548736787 357.09763745092783 25.490253793203657 268.023521027068 3.0
0.58904198203278091 0.72270060278093695 0.31942765324679046 350.68488850223423 130.4171465853421 250.42731366202318 4.0
0.51285660590703486 0.96889097226822973 0.65275467737350745 23.745542919457275 118.98401463018114 322.60963659419878 5.0
0.78608003485028433 0.83273743685098622 0.46591785719509976 124.52498788960992 100.66865249263579 43.350904777210218 6.0
0.65676045955005913 0.90612854270261067 0.46812684725311626 206.73481508655914 108.36640892186001 80.109515277983789 7.0
0.41091744799856139 0.019203430085754657 0.87577849258950335 294.38492822136715 146.40525644850072 307.47368257125362 8.0
0.2895339668620191 0.44890615451191845 0.98331278676555256 155.95129760119522 47.149690499466338 129.03566717283138 9.0
0.19961281156351873 0.52634383062850942 0.65188451822931848 147.12314868626314 111.70076966247582 118.18572187802707 10.0
0.86414247862963223 0.1358065510164656 0.66025345324864337 164.3847245485006 106.948282223783 169.81246394416348 11.0
0.22971651291623074 0.092972318577821886 0.29406405983067813 152.69170803150587 154.25570085621541 12.482717398044327 12.0
0.26338815658881415 0.34338560362947429 0.55845211616339796 34.576603888911734 112.1396081205236 231.97898012368159 13.0
0.75109304237913643 0.32426372309630619 0.24464858180476037 287.27773986438422 132.7748719439447 29.566044111233396 14.0
0.011464166371603362 0.038504815611266896 0.31848008962612995 3.6027692030412783 128.19004192002171 318.21386202740894 15.0
0.40531294455896061 0.89392258706810201 0.47360685251709117 224.94453046189483 91.073774858498993 174.6238603309032 16.0
0.53642882463725594 0.12961813440684475 0.33670742966203715 275.10050328051165 143.71902154901966 46.372591362351443 17.0
0.025264257063423813 0.86284946730733791 0.67853751997904233 286.09297442950589 84.366012495567063 168.12310601585438 18.0
0.46082042086486502 0.79920741984567956 0.84550103531963372 338.58981410067844 115.61172937509538 33.588172611417498 19.0
0.22570807057805362 0.074166418124772107 0.35703686595525042 123.22376691705952 84.092264279947017 358.5702863996658 20.0
0.05386086781200651 0.33174190751238741 0.22207351758975458 347.73707141532731 68.522081814108546 343.42676588519805 21.0
0.843158604433492 0.92955496315098074 0.64647123931005734 11.343815482295781 80.300931773797004 9.6393328996438079 22.0
0.38975306778625629 0.24157610260940071 0.71161594028191588 321.39703457206355 30.680985581522023 310.97284763119887 23.0
0.29080297238998321 0.7438587097696947 0.27827316089105131 318.66484094014749 129.93793511237541 136.82657482859585 24.0
0.39382389364070247 0.28978401907200979 0.25701142568390795 322.47065731551987 13.846167927307052 301.54027053054892 25.0
0.61050322346481545 0.13737535992809438 0.36661645869662263 352.54143971537871 57.8511858353625 133.84653788992898 26.0
0.79736663927764695 0.20513299822009629 0.79699332479250651 290.58637400802854 44.449209602954802 275.77563923277597 27.0
0.75235587126626513 0.11041486201059918 0.8131872750127791 70.389885527768058 106.61781772242031 249.0896396040977 28.0
0.47139010668774128 0.12192484253468709 0.21955576044612418 82.523861430871293 130.07642048077489 161.94830004765717 29.0
0.58577411200822327 0.55808726366080907 0.68861538513192688 4.5456602316904782 68.430488072013802 279.06105056042912 30.0
0.078221348390348527 0.38485150106633381 0.70002412594863284 44.840105036355524 52.915732353957182 321.10892793267385 31.0
0.67648574989589816 0.36189363050547918 0.1744438641736718 56.290857666353922 79.852422734452261 218.87802771695559 32.0
0.66993786328789628 0.24839196429109262 0.22913111586511459 90.545592617209479 111.73679898243722 50.777738624812869 33.0
0.97253038612350284 0.5008359837170796 0.22908814679929382 258.2784447839781 81.324197699117292 308.75839223966972 34.0
0.57267221923324418 0.57812183688041852 0.27747089968489891 44.241276881211661 104.39672542923724 263.41942696808212 35.0
0.20684173793886379 0.43993013267805814 0.65735383309297513 343.60408990114365 51.644327943351122 302.98734797140071 36.0
0.74510273339709676 0.73117975286639059 0.88155543772031653 318.38483613589898 93.903589849536274 302.06468871599935 37.0
0.96140945332061889 0.16540946028864878 0.40824265860818898 97.086714635901274 130.50888029759304 221.78895191070089 38.0
0.76663076605317781 0.85911002545479809 0.11281299879667539 163.06393615448818 43.363447677950042 338.05013375241901 39.0
0.41268673658765898 0.24787882796675886 0.57686480644197569 200.12920794363012 45.222523931505947 280.23271113977307 40.0
0.77256877568016891 0.88174830744168597 0.85149237688892054 116.81358850313981 71.413890894473454 115.54962789790765 41.0
0.26725724981852333 0.2962688497890511 0.89524301333622525 254.14781916777747 83.176346219908254 33.979304092964192 42.0
0.58047025880020098 0.57494408407976194 0.61595960318628096 334.70268656247265 42.480438737564974 177.92796756121371 43.0
0.52102440567302477 0.7145666401672387 0.21858506378351775 178.43052543384653 153.21174542887405 324.42119289220273 44.0
0.77321583279723483 0.96647383074249249 0.5062943967878929 230.42797261926012 99.507340620849902 169.75007570059978 45.0
0.3364367026326 0.45790436703027437 0.27197669375839439 218.70321774431869 60.819721511735267 217.80859716828817 46.0
0.41823530342173082 0.077759964416919514 0.66113722050248613 189.26108507623661 50.425749120256064 78.019878648192815 47.0
0.8754300454839713 0.094969845269609401 0.42632522145904467 250.899467172654 33.14582034295529 150.05888748377424 48.0
0.1950290416819265 0.59474264558516909 0.93298429220138601 232.236367110732 47.258083025548189 34.83912199551915 49.0
0.91993054481220637 0.48586729788450678 0.10933899155043697 246.05124283375034 131.539860458254 249.58739755697601 50.0

27
src/CMakeLists.txt
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
@ -15,15 +15,16 @@ add_dependencies(SYSTEM_ROUTINES C_ROUTINES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:SYSTEM_ROUTINES>)
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)
add_dependencies(DAMASK_INTERFACE PREC SYSTEM_ROUTINES)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>)
add_library(IO OBJECT "IO.f90")
@ -38,6 +39,10 @@ add_library(DEBUG OBJECT "debug.f90")
add_dependencies(DEBUG NUMERICS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DEBUG>)
add_library(CONFIG OBJECT "config.f90")
add_dependencies(CONFIG DEBUG)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:CONFIG>)
add_library(FEsolving OBJECT "FEsolving.f90")
add_dependencies(FEsolving DEBUG)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEsolving>)
@ -47,11 +52,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>)
@ -61,7 +66,7 @@ elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM")
endif()
add_library(MATERIAL OBJECT "material.f90")
add_dependencies(MATERIAL MESH)
add_dependencies(MATERIAL MESH CONFIG)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATERIAL>)
add_library(DAMASK_HELPERS OBJECT "lattice.f90")
@ -158,7 +163,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 +175,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)

13
src/CPFEM.f90
View File

@ -62,16 +62,18 @@ subroutine CPFEM_initAll(el,ip)
numerics_init
use debug, only: &
debug_init
use config, only: &
config_init
use FEsolving, only: &
FE_init
use math, only: &
math_init
use mesh, only: &
mesh_init
use lattice, only: &
lattice_init
use material, only: &
material_init
use lattice, only: &
lattice_init
use constitutive, only: &
constitutive_init
use crystallite, only: &
@ -93,6 +95,7 @@ subroutine CPFEM_initAll(el,ip)
call IO_init
call numerics_init
call debug_init
call config_init
call math_init
call FE_init
call mesh_init(ip, el) ! pass on coordinates to alter calcMode of first ip
@ -143,7 +146,8 @@ subroutine CPFEM_init
material_phase, &
homogState, &
phase_plasticity, &
plasticState, &
plasticState
use config, only: &
material_Nhomogenization
use crystallite, only: &
crystallite_F0, &
@ -310,7 +314,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
thermal_type, &
THERMAL_conduction_ID, &
phase_Nsources, &
material_homog, &
material_homog
use config, only: &
material_Nhomogenization
use crystallite, only: &
crystallite_partionedF,&

13
src/CPFEM2.f90
View File

@ -27,16 +27,18 @@ subroutine CPFEM_initAll(el,ip)
numerics_init
use debug, only: &
debug_init
use config, only: &
config_init
use FEsolving, only: &
FE_init
use math, only: &
math_init
use mesh, only: &
mesh_init
use lattice, only: &
lattice_init
use material, only: &
material_init
use lattice, only: &
lattice_init
use constitutive, only: &
constitutive_init
use crystallite, only: &
@ -64,6 +66,7 @@ subroutine CPFEM_initAll(el,ip)
#endif
call numerics_init
call debug_init
call config_init
call math_init
call FE_init
call mesh_init(ip, el) ! pass on coordinates to alter calcMode of first ip
@ -108,7 +111,8 @@ subroutine CPFEM_init
material_phase, &
homogState, &
phase_plasticity, &
plasticState, &
plasticState
use config, only: &
material_Nhomogenization
use crystallite, only: &
crystallite_F0, &
@ -228,7 +232,8 @@ subroutine CPFEM_age()
hydrogenfluxState, &
material_phase, &
phase_plasticity, &
phase_Nsources, &
phase_Nsources
use config, only: &
material_Nhomogenization
use crystallite, only: &
crystallite_partionedF,&

12
src/DAMASK_spectral.f90
View File

@ -355,8 +355,8 @@ program DAMASK_spectral
select case (loadCases(1)%ID(field))
case(FIELD_MECH_ID)
select case (spectral_solver)
case (DAMASK_spectral_SolverBasicPETSc_label)
call basicPETSc_init
case (DAMASK_spectral_SolverBasic_label)
call basic_init
case (DAMASK_spectral_SolverPolarisation_label)
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
@ -513,8 +513,8 @@ program DAMASK_spectral
select case(loadCases(currentLoadCase)%ID(field))
case(FIELD_MECH_ID)
select case (spectral_solver)
case (DAMASK_spectral_SolverBasicPETSc_label)
call BasicPETSc_forward (&
case (DAMASK_spectral_SolverBasic_label)
call Basic_forward (&
guess,timeinc,timeIncOld,remainingLoadCaseTime, &
deformation_BC = loadCases(currentLoadCase)%deformation, &
stress_BC = loadCases(currentLoadCase)%stress, &
@ -542,8 +542,8 @@ program DAMASK_spectral
select case(loadCases(currentLoadCase)%ID(field))
case(FIELD_MECH_ID)
select case (spectral_solver)
case (DAMASK_spectral_SolverBasicPETSc_label)
solres(field) = BasicPETSC_solution (&
case (DAMASK_spectral_SolverBasic_label)
solres(field) = Basic_solution (&
incInfo,timeinc,timeIncOld, &
stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rotation)

38
src/IO.f90
View File

@ -900,10 +900,10 @@ function IO_spotTagInPart(fileUnit,part,tag,Nsections)
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
foundNextPart: if (IO_getTag(line,'<','>') /= '') then
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
endif foundNextPart
if (IO_getTag(line,'[',']') /= '') section = section + 1_pInt ! found [section] identifier
if (section > 0_pInt) then
chunkPos = IO_stringPos(line)
@ -925,13 +925,10 @@ logical function IO_globalTagInPart(fileUnit,part,tag)
character(len=*),intent(in) :: part, & !< part in which tag is searched for
tag !< tag to search for
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: section
character(len=65536) :: line
IO_globalTagInPart = .false. ! assume to nowhere spot tag
section = 0_pInt
line =''
rewind(fileUnit)
@ -942,16 +939,20 @@ logical function IO_globalTagInPart(fileUnit,part,tag)
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
foundNextPart: if (IO_getTag(line,'<','>') /= '') then
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') section = section + 1_pInt ! found [section] identifier
if (section == 0_pInt) then
endif foundNextPart
foundFirstSection: if (IO_getTag(line,'[',']') /= '') then
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif foundFirstSection
chunkPos = IO_stringPos(line)
if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) & ! match
match: if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) then
IO_globalTagInPart = .true.
endif
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif match
enddo
end function IO_globalTagInPart
@ -981,6 +982,10 @@ pure function IO_stringPos(string)
if ( string(left:left) == '#' ) exit
IO_stringPos = [IO_stringPos,int(left, pInt), int(right, pInt)]
IO_stringPos(1) = IO_stringPos(1)+1_pInt
endOfString: if (right < left) then
IO_stringPos(IO_stringPos(1)*2+1) = len_trim(string)
exit
endif endOfString
enddo
end function IO_stringPos
@ -1545,6 +1550,17 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
case (136_pInt)
msg = 'zero entry on stiffness diagonal for transformed phase'
!--------------------------------------------------------------------------------------------------
! errors related to the parsing of material.config
case (140_pInt)
msg = 'key not found'
case (141_pInt)
msg = 'number of chunks in string differs'
case (142_pInt)
msg = 'empty list'
case (143_pInt)
msg = 'no value found for key'
!--------------------------------------------------------------------------------------------------
! material error messages and related messages in mesh
case (150_pInt)

1
src/commercialFEM_fileList.f90
View File

@ -6,6 +6,7 @@
#include "IO.f90"
#include "numerics.f90"
#include "debug.f90"
#include "config.f90"
#include "math.f90"
#include "FEsolving.f90"
#include "mesh.f90"

641
src/config.f90 Normal file
View File

@ -0,0 +1,641 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Reads in the material configuration from file
!> @details Reads the material configuration file, where solverJobName.materialConfig takes
!! precedence over material.config. Stores the raw strings and the positions of delimiters for the
!! parts 'homogenization', 'crystallite', 'phase', 'texture', and 'microstucture'
!--------------------------------------------------------------------------------------------------
module config
use prec, only: &
pReal, &
pInt
implicit none
private
type, private :: tPartitionedString
character(len=:), allocatable :: val
integer(pInt), dimension(:), allocatable :: pos
end type tPartitionedString
type, public :: tPartitionedStringList
type(tPartitionedString) :: string
type(tPartitionedStringList), pointer :: next => null()
contains
procedure :: add => add
procedure :: show => show
procedure :: keyExists => keyExists
procedure :: countKeys => countKeys
procedure :: getFloat => getFloat
procedure :: getInt => getInt
procedure :: getString => getString
procedure :: getFloats => getFloats
procedure :: getInts => getInts
procedure :: getStrings => getStrings
end type tPartitionedStringList
type(tPartitionedStringList), public :: emptyList
type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: &
phaseConfig, &
microstructureConfig, &
homogenizationConfig, &
textureConfig, &
crystalliteConfig
character(len=64), dimension(:), allocatable, public, protected :: &
phase_name, & !< name of each phase
homogenization_name, & !< name of each homogenization
crystallite_name, & !< name of each crystallite setting
microstructure_name, & !< name of each microstructure
texture_name !< name of each texture
! ToDo: make private, no one needs to know that
character(len=*), parameter, public :: &
MATERIAL_partHomogenization = 'homogenization', & !< keyword for homogenization part
MATERIAL_partCrystallite = 'crystallite', & !< keyword for crystallite part
MATERIAL_partPhase = 'phase', & !< keyword for phase part
MATERIAL_partMicrostructure = 'microstructure', & !< keyword for microstructure part
MATERIAL_partTexture = 'texture' !< keyword for texture part
! ToDo: Remove, use size(phaseConfig) etc
integer(pInt), public, protected :: &
material_Ntexture, & !< number of textures
material_Nphase, & !< number of phases
material_Nhomogenization, & !< number of homogenizations
material_Nmicrostructure, & !< number of microstructures
material_Ncrystallite !< number of crystallite settings
! ToDo: make private, no one needs to know that
character(len=*), parameter, public :: &
MATERIAL_configFile = 'material.config', & !< generic name for material configuration file
MATERIAL_localFileExt = 'materialConfig' !< extension of solver job name depending material configuration file
public :: config_init
contains
subroutine config_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_error, &
IO_open_file, &
IO_read, &
IO_lc, &
IO_open_jobFile_stat, &
IO_getTag, &
IO_timeStamp, &
IO_EOF
use debug, only: &
debug_level, &
debug_material, &
debug_levelBasic
implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt
integer(pInt) :: myDebug
character(len=65536) :: &
line, &
part
myDebug = debug_level(debug_material)
write(6,'(/,a)') ' <<<+- material init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ...open material.config file
rewind(fileUnit)
line = '' ! to have it initialized
do while (trim(line) /= IO_EOF)
part = IO_lc(IO_getTag(line,'<','>'))
select case (trim(part))
case (trim(material_partPhase))
call parseFile(line,phase_name,phaseConfig,FILEUNIT)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6)
case (trim(material_partMicrostructure))
call parseFile(line,microstructure_name,microstructureConfig,FILEUNIT)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed'; flush(6)
case (trim(material_partCrystallite))
call parseFile(line,crystallite_name,crystalliteConfig,FILEUNIT)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6)
case (trim(material_partHomogenization))
call parseFile(line,homogenization_name,homogenizationConfig,FILEUNIT)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed'; flush(6)
case (trim(material_partTexture))
call parseFile(line,texture_name,textureConfig,FILEUNIT)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6)
case default
line = IO_read(fileUnit)
end select
enddo
material_Nhomogenization = size(homogenizationConfig)
if (material_Nhomogenization < 1_pInt) call IO_error(160_pInt,ext_msg=material_partHomogenization)
material_Nmicrostructure = size(microstructureConfig)
if (material_Nmicrostructure < 1_pInt) call IO_error(160_pInt,ext_msg=material_partMicrostructure)
material_Ncrystallite = size(crystalliteConfig)
if (material_Ncrystallite < 1_pInt) call IO_error(160_pInt,ext_msg=material_partCrystallite)
material_Nphase = size(phaseConfig)
if (material_Nphase < 1_pInt) call IO_error(160_pInt,ext_msg=material_partPhase)
material_Ntexture = size(textureConfig)
if (material_Ntexture < 1_pInt) call IO_error(160_pInt,ext_msg=material_partTexture)
end subroutine config_init
!--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part in the material configuration file
!--------------------------------------------------------------------------------------------------
subroutine parseFile(line,&
sectionNames,part,fileUnit)
use IO, only: &
IO_read, &
IO_error, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringValue, &
IO_stringPos, &
IO_EOF
implicit none
integer(pInt), intent(in) :: fileUnit
character(len=*), dimension(:), allocatable, intent(inout) :: sectionNames
type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part
character(len=65536),intent(out) :: line
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: s
character(len=65536) :: devNull
character(len=64) :: tag
logical :: echo
echo = .false.
allocate(part(0))
s = 0_pInt
do while (trim(line) /= IO_EOF) ! read through sections of material part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
foundNextPart: if (IO_getTag(line,'<','>') /= '') then
devNull = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif foundNextPart
nextSection: if (IO_getTag(line,'[',']') /= '') then
s = s + 1_pInt
part = [part, emptyList]
tag = IO_getTag(line,'[',']')
GfortranBug86033: if (.not. allocated(sectionNames)) then
allocate(sectionNames(1),source=tag)
else GfortranBug86033
sectionNames = [sectionNames,tag]
endif GfortranBug86033
endif nextSection
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(trim(line),chunkPos,1_pInt)) ! extract key
inSection: if (s > 0_pInt) then
call part(s)%add(IO_lc(trim(line)))
else inSection
echo = (trim(tag) == '/echo/')
endif inSection
enddo
if (echo) then
do s = 1, size(sectionNames)
call part(s)%show()
end do
end if
end subroutine parseFile
!--------------------------------------------------------------------------------------------------
!> @brief add element
!> @details Adds a string together with the start/end position of chunks in this string. The new
!! element is added at the end of the list. Empty strings are not added. All strings are converted
!! to lower case
!--------------------------------------------------------------------------------------------------
subroutine add(this,string)
use IO, only: &
IO_isBlank, &
IO_lc, &
IO_stringPos
implicit none
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: string
type(tPartitionedStringList), pointer :: new, item
if (IO_isBlank(string)) return
allocate(new)
new%string%val = IO_lc (trim(string))
new%string%pos = IO_stringPos(trim(string))
item => this
do while (associated(item%next))
item => item%next
enddo
item%next => new
end subroutine add
!--------------------------------------------------------------------------------------------------
!> @brief prints all elements
!> @details Strings are printed in order of insertion (FIFO)
!--------------------------------------------------------------------------------------------------
subroutine show(this)
implicit none
class(tPartitionedStringList) :: this
type(tPartitionedStringList), pointer :: item
item => this%next
do while (associated(item))
write(6,'(a)') trim(item%string%val)
item => item%next
end do
end subroutine show
!--------------------------------------------------------------------------------------------------
!> @brief deallocates all elements of a given list
!> @details Strings are printed in order of insertion (FIFO)
!--------------------------------------------------------------------------------------------------
! subroutine free_all()
! implicit none
!
! type(node), pointer :: item
!
! do
! item => first
!
! if (associated(item) .eqv. .FALSE.) exit
!
! first => first%next
! deallocate(item)
! end do
! end subroutine free_all
!--------------------------------------------------------------------------------------------------
!> @brief reports wether a given key (string value at first position) exists in the list
!--------------------------------------------------------------------------------------------------
logical function keyExists(this,key)
use IO, only: &
IO_stringValue
implicit none
class(tPartitionedStringList), intent(in) :: this
character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item
keyExists = .false.
item => this%next
do while (associated(item) .and. .not. keyExists)
keyExists = trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)
item => item%next
end do
end function keyExists
!--------------------------------------------------------------------------------------------------
!> @brief count number of key appearances
!> @details traverses list and counts each occurrence of specified key
!--------------------------------------------------------------------------------------------------
integer(pInt) function countKeys(this,key)
use IO, only: &
IO_stringValue
implicit none
class(tPartitionedStringList), intent(in) :: this
character(len=*), intent(in) :: key
type(tPartitionedStringList), pointer :: item
countKeys = 0_pInt
item => this%next
do while (associated(item))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) &
countKeys = countKeys + 1_pInt
item => item%next
end do
end function countKeys
!--------------------------------------------------------------------------------------------------
!> @brief gets float value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given
!--------------------------------------------------------------------------------------------------
real(pReal) function getFloat(this,key,defaultVal)
use IO, only : &
IO_error, &
IO_stringValue, &
IO_FloatValue
implicit none
class(tPartitionedStringList), intent(in) :: this
character(len=*), intent(in) :: key
real(pReal), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
logical :: found
if (present(defaultVal)) getFloat = defaultVal
found = present(defaultVal)
item => this%next
do while (associated(item))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
getFloat = IO_FloatValue(item%string%val,item%string%pos,2)
endif
item => item%next
end do
if (.not. found) call IO_error(140_pInt,ext_msg=key)
end function getFloat
!--------------------------------------------------------------------------------------------------
!> @brief gets integer value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given
!--------------------------------------------------------------------------------------------------
integer(pInt) function getInt(this,key,defaultVal)
use IO, only: &
IO_error, &
IO_stringValue, &
IO_IntValue
implicit none
class(tPartitionedStringList), intent(in) :: this
character(len=*), intent(in) :: key
integer(pInt), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
logical :: found
if (present(defaultVal)) getInt = defaultVal
found = present(defaultVal)
item => this%next
do while (associated(item))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
getInt = IO_IntValue(item%string%val,item%string%pos,2)
endif
item => item%next
end do
if (.not. found) call IO_error(140_pInt,ext_msg=key)
end function getInt
!--------------------------------------------------------------------------------------------------
!> @brief gets string value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given. If raw is true, the the complete string is returned, otherwise
!! the individual chunks are returned
!--------------------------------------------------------------------------------------------------
character(len=65536) function getString(this,key,defaultVal,raw)
use IO, only: &
IO_error, &
IO_stringValue
implicit none
class(tPartitionedStringList), intent(in) :: this
character(len=*), intent(in) :: key
character(len=65536), intent(in), optional :: defaultVal
logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item
logical :: found, &
split
if (present(defaultVal)) getString = defaultVal
split = merge(.not. raw,.true.,present(raw))
found = present(defaultVal)
item => this%next
do while (associated(item))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
if (split) then
getString = IO_StringValue(item%string%val,item%string%pos,2)
else
getString = trim(item%string%val(item%string%pos(4):)) ! raw string starting a second chunk
endif
endif
item => item%next
end do
if (.not. found) call IO_error(140_pInt,ext_msg=key)
end function getString
!--------------------------------------------------------------------------------------------------
!> @brief gets array of float values of for a given key from a linked list
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getFloats(this,key,defaultVal)
use IO, only: &
IO_error, &
IO_stringValue, &
IO_FloatValue
implicit none
real(pReal), dimension(:), allocatable :: getFloats
class(tPartitionedStringList), intent(in) :: this
character(len=*), intent(in) :: key
real(pReal), dimension(:), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
integer(pInt) :: i
logical :: found, &
cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
found = .false.
allocate(getFloats(0))
item => this%next
do while (associated(item))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (.not. cumulative) then
deallocate(getFloats) ! use here rhs allocation with empty list
allocate(getFloats(0))
endif
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
do i = 2_pInt, item%string%pos(1)
getFloats = [getFloats,IO_FloatValue(item%string%val,item%string%pos,i)]
enddo
endif
item => item%next
end do
if (present(defaultVal) .and. .not. found) then
getFloats = defaultVal
found = .true.
endif
if (.not. found) call IO_error(140_pInt,ext_msg=key)
end function getFloats
!--------------------------------------------------------------------------------------------------
!> @brief gets array of integer values of for a given key from a linked list
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getInts(this,key,defaultVal)
use IO, only: &
IO_error, &
IO_stringValue, &
IO_IntValue
implicit none
integer(pInt), dimension(:), allocatable :: getInts
class(tPartitionedStringList), intent(in) :: this
character(len=*), intent(in) :: key
integer(pInt), dimension(:), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item
integer(pInt) :: i
logical :: found, &
cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
found = .false.
allocate(getInts(0))
item => this%next
do while (associated(item))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (.not. cumulative) then
deallocate(getInts) ! use here rhs allocation with empty list
allocate(getInts(0))
endif
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
do i = 2_pInt, item%string%pos(1)
getInts = [getInts,IO_IntValue(item%string%val,item%string%pos,i)]
enddo
endif
item => item%next
end do
if (present(defaultVal) .and. .not. found) then
getInts = defaultVal
found = .true.
endif
if (.not. found) call IO_error(140_pInt,ext_msg=key)
end function getInts
!--------------------------------------------------------------------------------------------------
!> @brief gets array of string values of for a given key from a linked list
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!! If raw is true, the the complete string is returned, otherwise the individual chunks are returned
!--------------------------------------------------------------------------------------------------
function getStrings(this,key,defaultVal,raw)
use IO, only: &
IO_error, &
IO_StringValue
implicit none
character(len=65536),dimension(:), allocatable :: getStrings
class(tPartitionedStringList), intent(in) :: this
character(len=*), intent(in) :: key
character(len=65536),dimension(:), intent(in), optional :: defaultVal
logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item
character(len=65536) :: str
integer(pInt) :: i
logical :: found, &
split, &
cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
split = merge(.not. raw,.true.,present(raw))
found = .false.
item => this%next
do while (associated(item))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true.
if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings)
if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
notAllocated: if (.not. allocated(getStrings)) then
if (split) then
str = IO_StringValue(item%string%val,item%string%pos,2_pInt)
allocate(getStrings(1),source=str)
do i=3_pInt,item%string%pos(1)
str = IO_StringValue(item%string%val,item%string%pos,i)
getStrings = [getStrings,str]
enddo
else
str = item%string%val(item%string%pos(4):)
getStrings = [str]
endif
else notAllocated
if (split) then
do i=2_pInt,item%string%pos(1)
str = IO_StringValue(item%string%val,item%string%pos,i)
getStrings = [getStrings,str]
enddo
else
getStrings = [getStrings,str]
endif
endif notAllocated
endif
item => item%next
end do
if (present(defaultVal) .and. .not. found) then
getStrings = defaultVal
found = .true.
endif
if (.not. found) call IO_error(140_pInt,ext_msg=key)
end function getStrings
end module config

32
src/constitutive.f90
View File

@ -59,12 +59,13 @@ subroutine constitutive_init()
IO_timeStamp
use mesh, only: &
FE_geomtype
use material, only: &
material_phase, &
use config, only: &
material_Nphase, &
material_localFileExt, &
material_configFile, &
phase_name, &
material_configFile
use material, only: &
material_phase, &
phase_plasticity, &
phase_plasticityInstance, &
phase_Nsources, &
@ -143,7 +144,6 @@ subroutine constitutive_init()
ins !< instance of plasticity/source
integer(pInt), dimension(:,:), pointer :: thisSize
integer(pInt), dimension(:) , pointer :: thisNoutput
character(len=64), dimension(:,:), pointer :: thisOutput
character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready
logical :: knownPlasticity, knownSource, nonlocalConstitutionPresent
@ -157,7 +157,7 @@ subroutine constitutive_init()
!--------------------------------------------------------------------------------------------------
! parse plasticities from config file
if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT)
@ -205,37 +205,30 @@ subroutine constitutive_init()
plasticityType: select case(phase_plasticity(p))
case (PLASTICITY_NONE_ID) plasticityType
outputName = PLASTICITY_NONE_label
thisNoutput => null()
thisOutput => null()
thisSize => null()
case (PLASTICITY_ISOTROPIC_ID) plasticityType
outputName = PLASTICITY_ISOTROPIC_label
thisNoutput => plastic_isotropic_Noutput
thisOutput => plastic_isotropic_output
thisSize => plastic_isotropic_sizePostResult
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
outputName = PLASTICITY_PHENOPOWERLAW_label
thisNoutput => plastic_phenopowerlaw_Noutput
thisOutput => plastic_phenopowerlaw_output
thisSize => plastic_phenopowerlaw_sizePostResult
case (PLASTICITY_KINEHARDENING_ID) plasticityType
outputName = PLASTICITY_KINEHARDENING_label
thisNoutput => plastic_kinehardening_Noutput
thisOutput => plastic_kinehardening_output
thisSize => plastic_kinehardening_sizePostResult
case (PLASTICITY_DISLOTWIN_ID) plasticityType
outputName = PLASTICITY_DISLOTWIN_label
thisNoutput => plastic_dislotwin_Noutput
thisOutput => plastic_dislotwin_output
thisSize => plastic_dislotwin_sizePostResult
case (PLASTICITY_DISLOUCLA_ID) plasticityType
outputName = PLASTICITY_DISLOUCLA_label
thisNoutput => plastic_disloucla_Noutput
thisOutput => plastic_disloucla_output
thisSize => plastic_disloucla_sizePostResult
case (PLASTICITY_NONLOCAL_ID) plasticityType
outputName = PLASTICITY_NONLOCAL_label
thisNoutput => plastic_nonlocal_Noutput
thisOutput => plastic_nonlocal_output
thisSize => plastic_nonlocal_sizePostResult
case default plasticityType
@ -246,7 +239,8 @@ subroutine constitutive_init()
write(FILEUNIT,'(a)') '(plasticity)'//char(9)//trim(outputName)
if (phase_plasticity(p) /= PLASTICITY_NONE_ID) then
OutputPlasticityLoop: do o = 1_pInt,thisNoutput(ins)
OutputPlasticityLoop: do o = 1_pInt,size(thisOutput(:,ins))
if(len(trim(thisOutput(o,ins))) > 0_pInt) &
write(FILEUNIT,'(a,i4)') trim(thisOutput(o,ins))//char(9),thisSize(o,ins)
enddo OutputPlasticityLoop
endif
@ -257,55 +251,46 @@ subroutine constitutive_init()
case (SOURCE_thermal_dissipation_ID) sourceType
ins = source_thermal_dissipation_instance(p)
outputName = SOURCE_thermal_dissipation_label
thisNoutput => source_thermal_dissipation_Noutput
thisOutput => source_thermal_dissipation_output
thisSize => source_thermal_dissipation_sizePostResult
case (SOURCE_thermal_externalheat_ID) sourceType
ins = source_thermal_externalheat_instance(p)
outputName = SOURCE_thermal_externalheat_label
thisNoutput => source_thermal_externalheat_Noutput
thisOutput => source_thermal_externalheat_output
thisSize => source_thermal_externalheat_sizePostResult
case (SOURCE_damage_isoBrittle_ID) sourceType
ins = source_damage_isoBrittle_instance(p)
outputName = SOURCE_damage_isoBrittle_label
thisNoutput => source_damage_isoBrittle_Noutput
thisOutput => source_damage_isoBrittle_output
thisSize => source_damage_isoBrittle_sizePostResult
case (SOURCE_damage_isoDuctile_ID) sourceType
ins = source_damage_isoDuctile_instance(p)
outputName = SOURCE_damage_isoDuctile_label
thisNoutput => source_damage_isoDuctile_Noutput
thisOutput => source_damage_isoDuctile_output
thisSize => source_damage_isoDuctile_sizePostResult
case (SOURCE_damage_anisoBrittle_ID) sourceType
ins = source_damage_anisoBrittle_instance(p)
outputName = SOURCE_damage_anisoBrittle_label
thisNoutput => source_damage_anisoBrittle_Noutput
thisOutput => source_damage_anisoBrittle_output
thisSize => source_damage_anisoBrittle_sizePostResult
case (SOURCE_damage_anisoDuctile_ID) sourceType
ins = source_damage_anisoDuctile_instance(p)
outputName = SOURCE_damage_anisoDuctile_label
thisNoutput => source_damage_anisoDuctile_Noutput
thisOutput => source_damage_anisoDuctile_output
thisSize => source_damage_anisoDuctile_sizePostResult
case (SOURCE_vacancy_phenoplasticity_ID) sourceType
ins = source_vacancy_phenoplasticity_instance(p)
outputName = SOURCE_vacancy_phenoplasticity_label
thisNoutput => source_vacancy_phenoplasticity_Noutput
thisOutput => source_vacancy_phenoplasticity_output
thisSize => source_vacancy_phenoplasticity_sizePostResult
case (SOURCE_vacancy_irradiation_ID) sourceType
ins = source_vacancy_irradiation_instance(p)
outputName = SOURCE_vacancy_irradiation_label
thisNoutput => source_vacancy_irradiation_Noutput
thisOutput => source_vacancy_irradiation_output
thisSize => source_vacancy_irradiation_sizePostResult
case (SOURCE_vacancy_thermalfluc_ID) sourceType
ins = source_vacancy_thermalfluc_instance(p)
outputName = SOURCE_vacancy_thermalfluc_label
thisNoutput => source_vacancy_thermalfluc_Noutput
thisOutput => source_vacancy_thermalfluc_output
thisSize => source_vacancy_thermalfluc_sizePostResult
case default sourceType
@ -313,7 +298,8 @@ subroutine constitutive_init()
end select sourceType
if (knownSource) then
write(FILEUNIT,'(a)') '(source)'//char(9)//trim(outputName)
OutputSourceLoop: do o = 1_pInt,thisNoutput(ins)
OutputSourceLoop: do o = 1_pInt,size(thisOutput(:,ins))
if(len(trim(thisOutput(o,ins))) > 0_pInt) &
write(FILEUNIT,'(a,i4)') trim(thisOutput(o,ins))//char(9),thisSize(o,ins)
enddo OutputSourceLoop
endif

215
src/crystallite.f90
View File

@ -155,7 +155,6 @@ subroutine crystallite_init
math_I3, &
math_EulerToR, &
math_inv33, &
math_transpose33, &
math_mul33xx33, &
math_mul33x33
use FEsolving, only: &
@ -167,28 +166,19 @@ subroutine crystallite_init
mesh_maxNips, &
mesh_maxNipNeighbors
use IO, only: &
IO_read, &
IO_timeStamp, &
IO_open_jobFile_stat, &
IO_open_file, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_write_jobFile, &
IO_error, &
IO_EOF
IO_error
use material
use config
use constitutive, only: &
constitutive_initialFi, &
constitutive_microstructure ! derived (shortcut) quantities of given state
implicit none
integer(pInt), parameter :: &
FILEUNIT = 200_pInt
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt), parameter :: FILEUNIT=434_pInt
integer(pInt) :: &
c, & !< counter in integration point component loop
i, & !< counter in integration point loop
@ -200,12 +190,11 @@ subroutine crystallite_init
eMax, & !< maximum number of elements
nMax, & !< maximum number of ip neighbors
myNcomponents, & !< number of components at current IP
section = 0_pInt, &
mySize
character(len=65536), dimension(:), allocatable :: str
character(len=65536) :: &
tag = '', &
line= ''
tag = ''
write(6,'(/,a)') ' <<<+- crystallite init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -277,85 +266,68 @@ subroutine crystallite_init
allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), &
material_Ncrystallite), source=0_pInt)
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ...open material.config file
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partCrystallite) ! wind forward to <crystallite>
line = IO_read(FILEUNIT)
enddo
do while (trim(line) /= IO_EOF) ! read through sections of crystallite part
line = IO_read(FILEUNIT)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(FILEUNIT, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt
o = 0_pInt ! reset output counter
cycle ! skip to next line
endif
if (section > 0_pInt) then
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
o = o + 1_pInt
crystallite_output(o,section) = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
outputName: select case(crystallite_output(o,section))
do c = 1_pInt, material_Ncrystallite
#if defined(__GFORTRAN__)
str = ['GfortranBug86277']
str = crystalliteConfig(c)%getStrings('(output)',defaultVal=str)
if (str(1) == 'GfortranBug86277') str = [character(len=65536)::]
#else
str = crystalliteConfig(c)%getStrings('(output)',defaultVal=[character(len=65536)::])
#endif
do o = 1_pInt, size(str)
crystallite_output(o,c) = str(o)
outputName: select case(str(o))
case ('phase') outputName
crystallite_outputID(o,section) = phase_ID
crystallite_outputID(o,c) = phase_ID
case ('texture') outputName
crystallite_outputID(o,section) = texture_ID
crystallite_outputID(o,c) = texture_ID
case ('volume') outputName
crystallite_outputID(o,section) = volume_ID
crystallite_outputID(o,c) = volume_ID
case ('grainrotationx') outputName
crystallite_outputID(o,section) = grainrotationx_ID
crystallite_outputID(o,c) = grainrotationx_ID
case ('grainrotationy') outputName
crystallite_outputID(o,section) = grainrotationy_ID
crystallite_outputID(o,c) = grainrotationy_ID
case ('grainrotationz') outputName
crystallite_outputID(o,section) = grainrotationx_ID
crystallite_outputID(o,c) = grainrotationx_ID
case ('orientation') outputName
crystallite_outputID(o,section) = orientation_ID
crystallite_outputID(o,c) = orientation_ID
case ('grainrotation') outputName
crystallite_outputID(o,section) = grainrotation_ID
crystallite_outputID(o,c) = grainrotation_ID
case ('eulerangles') outputName
crystallite_outputID(o,section) = eulerangles_ID
crystallite_outputID(o,c) = eulerangles_ID
case ('defgrad','f') outputName
crystallite_outputID(o,section) = defgrad_ID
crystallite_outputID(o,c) = defgrad_ID
case ('fe') outputName
crystallite_outputID(o,section) = fe_ID
crystallite_outputID(o,c) = fe_ID
case ('fp') outputName
crystallite_outputID(o,section) = fp_ID
crystallite_outputID(o,c) = fp_ID
case ('fi') outputName
crystallite_outputID(o,section) = fi_ID
crystallite_outputID(o,c) = fi_ID
case ('lp') outputName
crystallite_outputID(o,section) = lp_ID
crystallite_outputID(o,c) = lp_ID
case ('li') outputName
crystallite_outputID(o,section) = li_ID
crystallite_outputID(o,c) = li_ID
case ('e') outputName
crystallite_outputID(o,section) = e_ID
crystallite_outputID(o,c) = e_ID
case ('ee') outputName
crystallite_outputID(o,section) = ee_ID
crystallite_outputID(o,c) = ee_ID
case ('p','firstpiola','1stpiola') outputName
crystallite_outputID(o,section) = p_ID
crystallite_outputID(o,c) = p_ID
case ('s','tstar','secondpiola','2ndpiola') outputName
crystallite_outputID(o,section) = s_ID
crystallite_outputID(o,c) = s_ID
case ('elasmatrix') outputName
crystallite_outputID(o,section) = elasmatrix_ID
crystallite_outputID(o,c) = elasmatrix_ID
case ('neighboringip') outputName
crystallite_outputID(o,section) = neighboringip_ID
crystallite_outputID(o,c) = neighboringip_ID
case ('neighboringelement') outputName
crystallite_outputID(o,section) = neighboringelement_ID
crystallite_outputID(o,c) = neighboringelement_ID
case default outputName
call IO_error(105_pInt,ext_msg=IO_stringValue(line,chunkPos,2_pInt)//' (Crystallite)')
call IO_error(105_pInt,ext_msg=tag//' (Crystallite)')
end select outputName
end select
endif
enddo
enddo
close(FILEUNIT)
do r = 1_pInt,material_Ncrystallite
do o = 1_pInt,crystallite_Noutput(r)
@ -537,7 +509,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
use math, only: &
math_inv33, &
math_identity2nd, &
math_transpose33, &
math_mul33x33, &
math_mul66x6, &
math_Mandel6to33, &
@ -616,17 +587,17 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
write(6,'(/,a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> boundary values at el ip ipc ', &
debug_e,'(',mesh_element(1,debug_e), ')',debug_i, debug_g
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> F ', &
math_transpose33(crystallite_partionedF(1:3,1:3,debug_g,debug_i,debug_e))
transpose(crystallite_partionedF(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> F0 ', &
math_transpose33(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e))
transpose(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp0', &
math_transpose33(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e))
transpose(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fi0', &
math_transpose33(crystallite_partionedFi0(1:3,1:3,debug_g,debug_i,debug_e))
transpose(crystallite_partionedFi0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp0', &
math_transpose33(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
transpose(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Li0', &
math_transpose33(crystallite_partionedLi0(1:3,1:3,debug_g,debug_i,debug_e))
transpose(crystallite_partionedLi0(1:3,1:3,debug_g,debug_i,debug_e))
endif
!--------------------------------------------------------------------------------------------------
@ -1107,15 +1078,15 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
.or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip ipc ',e,i,c
write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', &
math_transpose33(crystallite_P(1:3,1:3,c,i,e))*1.0e-6_pReal
transpose(crystallite_P(1:3,1:3,c,i,e))*1.0e-6_pReal
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', &
math_transpose33(crystallite_Fp(1:3,1:3,c,i,e))
transpose(crystallite_Fp(1:3,1:3,c,i,e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fi', &
math_transpose33(crystallite_Fi(1:3,1:3,c,i,e))
transpose(crystallite_Fi(1:3,1:3,c,i,e))
write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST >> Lp', &
math_transpose33(crystallite_Lp(1:3,1:3,c,i,e))
transpose(crystallite_Lp(1:3,1:3,c,i,e))
write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST >> Li', &
math_transpose33(crystallite_Li(1:3,1:3,c,i,e))
transpose(crystallite_Li(1:3,1:3,c,i,e))
flush(6)
endif
enddo
@ -1166,7 +1137,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
temp_33 = transpose(math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
crystallite_invFi(1:3,1:3,c,i,e)))
rhs_3333 = 0.0_pReal
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
@ -1208,12 +1179,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal
temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), &
math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e))))
transpose(crystallite_invFp(1:3,1:3,c,i,e))))
forall(p=1_pInt:3_pInt) &
crystallite_dPdF(p,1:3,p,1:3,c,i,e) = math_transpose33(temp_33)
crystallite_dPdF(p,1:3,p,1:3,c,i,e) = transpose(temp_33)
temp_33 = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)), &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))
transpose(crystallite_invFp(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33)
@ -1223,14 +1194,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(math_mul33x33(temp_33,dSdF(1:3,1:3,p,o)), &
math_transpose33(crystallite_invFp(1:3,1:3,c,i,e)))
transpose(crystallite_invFp(1:3,1:3,c,i,e)))
temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e)), &
math_Mandel6to33(crystallite_Tstar_v(1:6,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(temp_33,math_transpose33(dFpinvdF(1:3,1:3,p,o)))
math_mul33x33(temp_33,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo
enddo elementLooping6
@ -1272,8 +1243,9 @@ subroutine crystallite_integrateStateRK4()
plasticState, &
sourceState, &
phase_Nsources, &
material_Nphase, &
phaseAt, phasememberAt
use config, only: &
material_Nphase
use constitutive, only: &
constitutive_collectDotState, &
constitutive_microstructure
@ -3195,7 +3167,6 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33)
use math, only: &
math_mul33x33, &
math_inv33, &
math_transpose33, &
math_EulerToR
use material, only: &
material_EulerAngles
@ -3210,8 +3181,8 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33)
ipc ! grain index
T = math_mul33x33(math_EulerToR(material_EulerAngles(1:3,ipc,ip,el)), &
math_transpose33(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el))))
crystallite_push33ToRef = math_mul33x33(math_transpose33(T),math_mul33x33(tensor33,T))
transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el))))
crystallite_push33ToRef = math_mul33x33(transpose(T),math_mul33x33(tensor33,T))
end function crystallite_push33ToRef
@ -3260,7 +3231,6 @@ logical function crystallite_integrateStress(&
math_mul3333xx3333, &
math_mul66x6, &
math_mul99x99, &
math_transpose33, &
math_inv33, &
math_invert, &
math_det33, &
@ -3386,7 +3356,7 @@ logical function crystallite_integrateStress(&
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip ipc ',&
el,'(',mesh_element(1,el),')',ip,ipc
if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0_pInt) &
write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',math_transpose33(Fp_current(1:3,1:3))
write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',transpose(Fp_current(1:3,1:3))
endif
#endif
return
@ -3402,7 +3372,7 @@ logical function crystallite_integrateStress(&
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip ipc ',&
el,'(',mesh_element(1,el),')',ip,ipc
if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0_pInt) &
write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',math_transpose33(Fi_current(1:3,1:3))
write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',transpose(Fi_current(1:3,1:3))
endif
#endif
return
@ -3465,9 +3435,9 @@ logical function crystallite_integrateStress(&
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i3,/)') '<< CRYST >> stress iteration ', NiterationStressLp
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Fi', math_transpose33(Fi_new)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Fe', math_transpose33(Fe)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Lpguess', transpose(Lpguess)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Fi', transpose(Fi_new)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Fe', transpose(Fe)
write(6,'(a,/,6(e20.10,1x))') '<< CRYST >> Tstar', Tstar_v
endif
#endif
@ -3488,7 +3458,7 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Lp_constitutive', transpose(Lp_constitutive)
endif
#endif
@ -3534,7 +3504,7 @@ logical function crystallite_integrateStress(&
if (mod(jacoCounterLp, iJacoLpresiduum) == 0_pInt) then
dFe_dLp3333 = 0.0_pReal
forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) &
dFe_dLp3333(o,1:3,p,1:3) = A(o,p)*math_transpose33(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
dFe_dLp3333(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
dFe_dLp3333 = - dt * dFe_dLp3333
dRLp_dLp = math_identity2nd(9_pInt) &
- math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333))
@ -3564,10 +3534,10 @@ logical function crystallite_integrateStress(&
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(math_Plain3333to99(dFe_dLp3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFe_constitutive',transpose(math_Plain3333to99(dT_dFe3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(math_Plain3333to99(dLp_dT3333))
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',math_transpose33(A)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',math_transpose33(B)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive',math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess',math_transpose33(Lpguess)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',transpose(A)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',transpose(B)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive',transpose(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess',transpose(Lpguess)
endif
endif
#endif
@ -3597,8 +3567,8 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive', math_transpose33(Li_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess', math_transpose33(Liguess)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive', transpose(Li_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess', transpose(Liguess)
endif
#endif
!* update current residuum and check for convergence of loop
@ -3653,8 +3623,8 @@ logical function crystallite_integrateStress(&
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333))
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',math_transpose33(Li_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',math_transpose33(Liguess)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',transpose(Li_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',transpose(Liguess)
endif
endif
#endif
@ -3688,7 +3658,7 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) &
write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new)
write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',transpose(invFp_new)
endif
#endif
return
@ -3699,7 +3669,7 @@ logical function crystallite_integrateStress(&
crystallite_P(1:3,1:3,ipc,ip,el) = math_mul33x33(math_mul33x33(Fg_new,invFp_new), &
math_mul33x33(math_Mandel6to33(Tstar_v), &
math_transpose33(invFp_new)))
transpose(invFp_new)))
!* store local values in global variables
@ -3719,13 +3689,13 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,ipc,ip,el))*1.0e-6_pReal
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',transpose(crystallite_P(1:3,1:3,ipc,ip,el))*1.0e-6_pReal
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Cauchy / MPa', &
math_mul33x33(crystallite_P(1:3,1:3,ipc,ip,el), math_transpose33(Fg_new)) * 1.0e-6_pReal / math_det33(Fg_new)
math_mul33x33(crystallite_P(1:3,1:3,ipc,ip,el), transpose(Fg_new)) * 1.0e-6_pReal / math_det33(Fg_new)
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fe Lp Fe^-1', &
math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,ipc,ip,el), math_inv33(Fe_new)))) ! transpose to get correct print out order
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,ipc,ip,el))
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fi',math_transpose33(crystallite_Fi(1:3,1:3,ipc,ip,el))
transpose(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,ipc,ip,el), math_inv33(Fe_new)))) ! transpose to get correct print out order
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp',transpose(crystallite_Fp(1:3,1:3,ipc,ip,el))
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fi',transpose(crystallite_Fi(1:3,1:3,ipc,ip,el))
endif
#endif
@ -3842,7 +3812,6 @@ function crystallite_postResults(ipc, ip, el)
math_qToEuler, &
math_qToEulerAxisAngle, &
math_mul33x33, &
math_transpose33, &
math_det33, &
math_I3, &
inDeg, &
@ -3945,41 +3914,41 @@ function crystallite_postResults(ipc, ip, el)
case (defgrad_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_partionedF(1:3,1:3,ipc,ip,el)),[mySize])
reshape(transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)),[mySize])
case (e_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( &
math_transpose33(crystallite_partionedF(1:3,1:3,ipc,ip,el)), &
transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)), &
crystallite_partionedF(1:3,1:3,ipc,ip,el)) - math_I3),[mySize])
case (fe_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_Fe(1:3,1:3,ipc,ip,el)),[mySize])
reshape(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)),[mySize])
case (ee_ID)
Ee = 0.5_pReal *(math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,ipc,ip,el)), &
Ee = 0.5_pReal *(math_mul33x33(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)), &
crystallite_Fe(1:3,1:3,ipc,ip,el)) - math_I3)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(Ee,[mySize])
case (fp_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_Fp(1:3,1:3,ipc,ip,el)),[mySize])
reshape(transpose(crystallite_Fp(1:3,1:3,ipc,ip,el)),[mySize])
case (fi_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_Fi(1:3,1:3,ipc,ip,el)),[mySize])
reshape(transpose(crystallite_Fi(1:3,1:3,ipc,ip,el)),[mySize])
case (lp_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_Lp(1:3,1:3,ipc,ip,el)),[mySize])
reshape(transpose(crystallite_Lp(1:3,1:3,ipc,ip,el)),[mySize])
case (li_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_Li(1:3,1:3,ipc,ip,el)),[mySize])
reshape(transpose(crystallite_Li(1:3,1:3,ipc,ip,el)),[mySize])
case (p_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(math_transpose33(crystallite_P(1:3,1:3,ipc,ip,el)),[mySize])
reshape(transpose(crystallite_P(1:3,1:3,ipc,ip,el)),[mySize])
case (s_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &

3
src/damage_local.f90
View File

@ -70,7 +70,8 @@ subroutine damage_local_init(fileUnit)
damageState, &
damageMapping, &
damage, &
damage_initialPhi, &
damage_initialPhi
use config, only: &
material_partHomogenization
implicit none

1
src/damage_none.f90
View File

@ -26,6 +26,7 @@ subroutine damage_none_init()
use IO, only: &
IO_timeStamp
use material
use config
implicit none
integer(pInt) :: &

3
src/damage_nonlocal.f90
View File

@ -75,7 +75,8 @@ subroutine damage_nonlocal_init(fileUnit)
damageState, &
damageMapping, &
damage, &
damage_initialPhi, &
damage_initialPhi
use config, only: &
material_partHomogenization
implicit none

5
src/homogenization.f90
View File

@ -101,6 +101,7 @@ subroutine homogenization_init
crystallite_maxSizePostResults
#endif
use material
use config
use homogenization_none
use homogenization_isostrain
use homogenization_RGC
@ -443,11 +444,9 @@ subroutine homogenization_init
allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems))
#endif
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
#ifdef TODO
@ -475,7 +474,7 @@ subroutine homogenization_init
flush(6)
if (debug_g < 1 .or. debug_g > homogenization_Ngrains(mesh_element(3,debug_e))) &
call IO_error(602_pInt,ext_msg='component (grain)', el=debug_e, g=debug_g)
call IO_error(602_pInt,ext_msg='constituent', el=debug_e, g=debug_g)
end subroutine homogenization_init

1
src/homogenization_RGC.f90
View File

@ -100,6 +100,7 @@ subroutine homogenization_RGC_init(fileUnit)
FE_geomtype
use IO
use material
use config
implicit none
integer(pInt), intent(in) :: fileUnit !< file pointer to material configuration

1
src/homogenization_isostrain.f90
View File

@ -62,6 +62,7 @@ subroutine homogenization_isostrain_init(fileUnit)
debug_levelBasic
use IO
use material
use config
implicit none
integer(pInt), intent(in) :: fileUnit

1
src/homogenization_none.f90
View File

@ -29,6 +29,7 @@ subroutine homogenization_none_init()
use IO, only: &
IO_timeStamp
use material
use config
implicit none
integer(pInt) :: &

3
src/hydrogenflux_cahnhilliard.f90
View File

@ -81,7 +81,8 @@ subroutine hydrogenflux_cahnhilliard_init(fileUnit)
hydrogenfluxMapping, &
hydrogenConc, &
hydrogenConcRate, &
hydrogenflux_initialCh, &
hydrogenflux_initialCh
use config, only: &
material_partHomogenization, &
material_partPhase

1
src/hydrogenflux_isoconc.f90
View File

@ -27,6 +27,7 @@ subroutine hydrogenflux_isoconc_init()
use IO, only: &
IO_timeStamp
use material
use config
implicit none
integer(pInt) :: &

3
src/kinematics_cleavage_opening.f90
View File

@ -78,7 +78,8 @@ subroutine kinematics_cleavage_opening_init(fileUnit)
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_cleavage_opening_label, &
KINEMATICS_cleavage_opening_ID, &
KINEMATICS_cleavage_opening_ID
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use lattice, only: &

3
src/kinematics_hydrogen_strain.f90
View File

@ -68,7 +68,8 @@ subroutine kinematics_hydrogen_strain_init(fileUnit)
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_hydrogen_strain_label, &
KINEMATICS_hydrogen_strain_ID, &
KINEMATICS_hydrogen_strain_ID
use config, only: &
material_Nphase, &
MATERIAL_partPhase

3
src/kinematics_slipplane_opening.f90
View File

@ -78,7 +78,8 @@ subroutine kinematics_slipplane_opening_init(fileUnit)
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_slipplane_opening_label, &
KINEMATICS_slipplane_opening_ID, &
KINEMATICS_slipplane_opening_ID
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use lattice, only: &

3
src/kinematics_thermal_expansion.f90
View File

@ -68,7 +68,8 @@ subroutine kinematics_thermal_expansion_init(fileUnit)
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_thermal_expansion_label, &
KINEMATICS_thermal_expansion_ID, &
KINEMATICS_thermal_expansion_ID
use config, only: &
material_Nphase, &
MATERIAL_partPhase

3
src/kinematics_vacancy_strain.f90
View File

@ -68,7 +68,8 @@ subroutine kinematics_vacancy_strain_init(fileUnit)
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_vacancy_strain_label, &
KINEMATICS_vacancy_strain_ID, &
KINEMATICS_vacancy_strain_ID
use config, only: &
material_Nphase, &
MATERIAL_partPhase

2
src/lattice.f90
View File

@ -1263,7 +1263,7 @@ subroutine lattice_init
IO_stringPos, &
IO_stringValue, &
IO_floatValue
use material, only: &
use config, only: &
material_configfile, &
material_localFileExt, &
material_partPhase

973
src/material.f90
View File

File diff suppressed because it is too large Load Diff

134
src/math.f90
View File

@ -223,7 +223,6 @@ end subroutine math_init
!> @brief check correctness of (some) math functions
!--------------------------------------------------------------------------------------------------
subroutine math_check
use prec, only: tol_math_check
use IO, only: IO_error
@ -1821,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))]
@ -2633,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

3
src/plastic_disloUCLA.f90
View File

@ -152,7 +152,8 @@ subroutine plastic_disloUCLA_init(fileUnit)
PLASTICITY_DISLOUCLA_label, &
PLASTICITY_DISLOUCLA_ID, &
material_phase, &
plasticState, &
plasticState
use config, only: &
MATERIAL_partPhase
use lattice
use numerics,only: &

3
src/plastic_dislotwin.f90
View File

@ -251,7 +251,8 @@ subroutine plastic_dislotwin_init(fileUnit)
PLASTICITY_DISLOTWIN_label, &
PLASTICITY_DISLOTWIN_ID, &
material_phase, &
plasticState, &
plasticState
use config, only: &
MATERIAL_partPhase
use lattice
use numerics,only: &

337
src/plastic_isotropic.f90
View File

@ -13,15 +13,10 @@ module plastic_isotropic
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_isotropic_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: &
plastic_isotropic_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_isotropic_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
plastic_isotropic_Noutput !< number of outputs per instance
@ -40,17 +35,17 @@ module plastic_isotropic
gdot0, &
n, &
h0, &
h0_slopeLnRate = 0.0_pReal, &
h0_slopeLnRate, &
tausat, &
a, &
aTolFlowstress = 1.0_pReal, &
aTolShear = 1.0e-6_pReal, &
tausat_SinhFitA= 0.0_pReal, &
tausat_SinhFitB= 0.0_pReal, &
tausat_SinhFitC= 0.0_pReal, &
tausat_SinhFitD= 0.0_pReal
aTolFlowstress, &
aTolShear, &
tausat_SinhFitA, &
tausat_SinhFitB, &
tausat_SinhFitC, &
tausat_SinhFitD
logical :: &
dilatation = .false.
dilatation
end type
type(tParameters), dimension(:), allocatable, target, private :: param !< containers of constitutive parameters (len Ninstance)
@ -79,12 +74,13 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_init(fileUnit)
subroutine plastic_isotropic_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO
use debug, only: &
debug_level, &
debug_constitutive, &
@ -94,17 +90,6 @@ subroutine plastic_isotropic_init(fileUnit)
use math, only: &
math_Mandel3333to66, &
math_Voigt66to3333
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_plasticity, &
phase_plasticityInstance, &
@ -112,17 +97,17 @@ subroutine plastic_isotropic_init(fileUnit)
PLASTICITY_ISOTROPIC_label, &
PLASTICITY_ISOTROPIC_ID, &
material_phase, &
plasticState, &
MATERIAL_partPhase
plasticState
use config, only: &
MATERIAL_partPhase, &
phaseConfig
use lattice
implicit none
integer(pInt), intent(in) :: fileUnit
type(tParameters), pointer :: p
type(tParameters), pointer :: prm
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: &
o, &
phase, &
@ -133,174 +118,103 @@ subroutine plastic_isotropic_init(fileUnit)
sizeState, &
sizeDeltaState
character(len=65536) :: &
tag = '', &
line = '', &
extmsg = ''
character(len=64) :: &
outputtag = ''
integer(pInt) :: NipcMyPhase
integer(pInt) :: NipcMyPhase,i
character(len=65536), dimension(:), allocatable :: outputs
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(phase_plasticity == PLASTICITY_ISOTROPIC_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(plastic_isotropic_sizePostResults(maxNinstance), source=0_pInt)
! public variables
allocate(plastic_isotropic_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
allocate(plastic_isotropic_output(maxval(phase_Noutput), maxNinstance))
plastic_isotropic_output = ''
allocate(plastic_isotropic_Noutput(maxNinstance), source=0_pInt)
! inernal variable
allocate(param(maxNinstance)) ! one container of parameters per instance
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
phase = phase + 1_pInt ! advance section counter
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
select case(tag)
case ('(output)')
outputtag = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
select case(outputtag)
case ('flowstress')
plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
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
plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag
p%outputID = [p%outputID,strainrate_ID]
end select
case ('/dilatation/')
p%dilatation = .true.
case ('tau0')
p%tau0 = IO_floatValue(line,chunkPos,2_pInt)
case ('gdot0')
p%gdot0 = IO_floatValue(line,chunkPos,2_pInt)
case ('n')
p%n = IO_floatValue(line,chunkPos,2_pInt)
case ('h0')
p%h0 = IO_floatValue(line,chunkPos,2_pInt)
case ('h0_slope','slopelnrate')
p%h0_slopeLnRate = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat')
p%tausat = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfita')
p%tausat_SinhFitA = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitb')
p%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitc')
p%tausat_SinhFitC = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitd')
p%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt)
case ('a', 'w0')
p%a = IO_floatValue(line,chunkPos,2_pInt)
case ('taylorfactor')
p%fTaylor = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_flowstress')
p%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_shear')
p%aTolShear = IO_floatValue(line,chunkPos,2_pInt)
case default
end select
endif; endif
enddo parsingFile
allocate(state(maxNinstance)) ! internal state aliases
allocate(dotState(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)
do phase = 1_pInt, size(phase_plasticityInstance)
if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then
instance = phase_plasticityInstance(phase)
p => param(instance)
extmsg = ''
prm => param(instance) ! shorthand pointer to parameter object of my constitutive law
prm%tau0 = phaseConfig(phase)%getFloat('tau0')
prm%tausat = phaseConfig(phase)%getFloat('tausat')
prm%gdot0 = phaseConfig(phase)%getFloat('gdot0')
prm%n = phaseConfig(phase)%getFloat('n')
prm%h0 = phaseConfig(phase)%getFloat('h0')
prm%fTaylor = phaseConfig(phase)%getFloat('m')
prm%h0_slopeLnRate = phaseConfig(phase)%getFloat('h0_slopelnrate', defaultVal=0.0_pReal)
prm%tausat_SinhFitA = phaseConfig(phase)%getFloat('tausat_sinhfita',defaultVal=0.0_pReal)
prm%tausat_SinhFitB = phaseConfig(phase)%getFloat('tausat_sinhfitb',defaultVal=0.0_pReal)
prm%tausat_SinhFitC = phaseConfig(phase)%getFloat('tausat_sinhfitc',defaultVal=0.0_pReal)
prm%tausat_SinhFitD = phaseConfig(phase)%getFloat('tausat_sinhfitd',defaultVal=0.0_pReal)
prm%a = phaseConfig(phase)%getFloat('a')
prm%aTolFlowStress = phaseConfig(phase)%getFloat('atol_flowstress',defaultVal=1.0_pReal)
prm%aTolShear = phaseConfig(phase)%getFloat('atol_shear',defaultVal=1.0e-6_pReal)
prm%dilatation = phaseConfig(phase)%keyExists('/dilatation/')
#if defined(__GFORTRAN__)
outputs = ['GfortranBug86277']
outputs = phaseConfig(phase)%getStrings('(output)',defaultVal=outputs)
if (outputs(1) == 'GfortranBug86277') outputs = [character(len=65536)::]
#else
outputs = phaseConfig(phase)%getStrings('(output)',defaultVal=[character(len=65536)::])
#endif
allocate(prm%outputID(0))
do i=1_pInt, size(outputs)
select case(outputs(i))
case ('flowstress')
plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputs(i)
plasticState(phase)%sizePostResults = plasticState(phase)%sizePostResults + 1_pInt
plastic_isotropic_sizePostResult(i,instance) = 1_pInt
prm%outputID = [prm%outputID,flowstress_ID]
case ('strainrate')
plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputs(i)
plasticState(phase)%sizePostResults = &
plasticState(phase)%sizePostResults + 1_pInt
plastic_isotropic_sizePostResult(i,instance) = 1_pInt
prm%outputID = [prm%outputID,strainrate_ID]
end select
enddo
!--------------------------------------------------------------------------------------------------
! sanity checks
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)
endif
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
select case(p%outputID(o))
case(flowstress_ID,strainrate_ID)
mySize = 1_pInt
case default
end select
outputFound: if (mySize > 0_pInt) then
plastic_isotropic_sizePostResult(o,instance) = mySize
plastic_isotropic_sizePostResults(instance) = &
plastic_isotropic_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
extmsg = ''
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//"'aTolShear' "
if (prm%tau0 < 0.0_pReal) extmsg = trim(extmsg)//"'tau0' "
if (prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//"'gdot0' "
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//"'n' "
if (prm%tausat <= prm%tau0) extmsg = trim(extmsg)//"'tausat' "
if (prm%a <= 0.0_pReal) extmsg = trim(extmsg)//"'a' "
if (prm%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//"'m' "
if (prm%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//"'atol_flowstress' "
if (extmsg /= '') call IO_error(211_pInt,ip=instance,&
ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')')
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc)
sizeDotState = size(["flowstress ","accumulated_shear"])
sizeDeltaState = 0_pInt ! no sudden jumps in state
sizeState = sizeDotState + sizeDeltaState
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%sizePostResults = plastic_isotropic_sizePostResults(instance)
plasticState(phase)%nSlip = 1
plasticState(phase)%nTwin = 0
plasticState(phase)%nTrans= 0
allocate(plasticState(phase)%aTolState ( sizeState))
allocate(plasticState(phase)%state0 ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%partionedState0 ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%subState0 ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%state ( sizeState,NipcMyPhase),source=0.0_pReal)
@ -320,22 +234,23 @@ subroutine plastic_isotropic_init(fileUnit)
state(instance)%flowstress => plasticState(phase)%state (1,1:NipcMyPhase)
dotState(instance)%flowstress => plasticState(phase)%dotState (1,1:NipcMyPhase)
plasticState(phase)%state0(1,1:NipcMyPhase) = p%tau0
plasticState(phase)%aTolState(1) = p%aTolFlowstress
plasticState(phase)%state0(1,1:NipcMyPhase) = prm%tau0
plasticState(phase)%aTolState(1) = prm%aTolFlowstress
state(instance)%accumulatedShear => plasticState(phase)%state (2,1:NipcMyPhase)
dotState(instance)%accumulatedShear => plasticState(phase)%dotState (2,1:NipcMyPhase)
plasticState(phase)%state0 (2,1:NipcMyPhase) = 0.0_pReal
plasticState(phase)%aTolState(2) = p%aTolShear
plasticState(phase)%aTolState(2) = prm%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
endif
enddo
end subroutine plastic_isotropic_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
@ -354,8 +269,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
math_Mandel6to33, &
math_Plain3333to99, &
math_deviatoric33, &
math_mul33xx33, &
math_transpose33
math_mul33xx33
use material, only: &
phasememberAt, &
material_phase, &
@ -374,7 +288,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
ip, & !< integration point
el !< element
type(tParameters), pointer :: p
type(tParameters), pointer :: prm
real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
@ -390,7 +304,7 @@ 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)
prm => 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)
@ -400,31 +314,31 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
Lp = 0.0_pReal
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = p%gdot0 &
* ( sqrt(1.5_pReal) * norm_Tstar_dev / p%fTaylor / state(instance)%flowstress(of) ) &
**p%n
gamma_dot = prm%gdot0 &
* ( sqrt(1.5_pReal) * norm_Tstar_dev / prm%fTaylor / state(instance)%flowstress(of) ) &
**prm%n
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/p%fTaylor
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/prm%fTaylor
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CONST isotropic >> at el ip g ',el,ip,ipc
write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CONST isotropic >> Tstar (dev) / MPa', &
math_transpose33(Tstar_dev_33(1:3,1:3))*1.0e-6_pReal
transpose(Tstar_dev_33(1:3,1:3))*1.0e-6_pReal
write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> norm Tstar / MPa', norm_Tstar_dev*1.0e-6_pReal
write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', gamma_dot
end if
!--------------------------------------------------------------------------------------------------
! 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) = (p%n-1.0_pReal) * &
dLp_dTstar_3333(k,l,m,n) = (prm%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 / p%fTaylor * &
dLp_dTstar99 = math_Plain3333to99(gamma_dot / prm%fTaylor * &
dLp_dTstar_3333 / norm_Tstar_dev)
end if
end subroutine plastic_isotropic_LpAndItsTangent
@ -456,7 +370,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
ip, & !< integration point
el !< element
type(tParameters), pointer :: p
type(tParameters), pointer :: prm
real(pReal), dimension(3,3) :: &
Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
@ -470,28 +384,28 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
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)
prm => 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 (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
if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero
gamma_dot = prm%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / prm%fTaylor / state(instance)%flowstress(of) ) &
**prm%n
Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/p%fTaylor
Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/prm%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) = (p%n-1.0_pReal) * &
dLi_dTstar_3333(k,l,m,n) = (prm%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 / p%fTaylor * &
dLi_dTstar_3333 = gamma_dot / prm%fTaylor * &
dLi_dTstar_3333 / norm_Tstar_sph
else
Li = 0.0_pReal
@ -520,7 +434,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
type(tParameters), pointer :: prm
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
@ -534,11 +448,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)
prm => param(instance)
!--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress
if (p%dilatation) then
if (prm%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
@ -547,26 +461,26 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
end if
!--------------------------------------------------------------------------------------------------
! strain rate
gamma_dot = p%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
gamma_dot = prm%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
/ &!-----------------------------------------------------------------------------------
(p%fTaylor*state(instance)%flowstress(of) ))**p%n
(prm%fTaylor*state(instance)%flowstress(of) ))**prm%n
!--------------------------------------------------------------------------------------------------
! hardening coefficient
if (abs(gamma_dot) > 1e-12_pReal) then
if (dEq0(p%tausat_SinhFitA)) then
saturation = p%tausat
if (dEq0(prm%tausat_SinhFitA)) then
saturation = prm%tausat
else
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) &
saturation = prm%tausat &
+ asinh( (gamma_dot / prm%tausat_SinhFitA&
)**(1.0_pReal / prm%tausat_SinhFitD)&
)**(1.0_pReal / prm%tausat_SinhFitC) &
/ ( prm%tausat_SinhFitB &
* (gamma_dot / prm%gdot0)**(1.0_pReal / prm%n) &
)
endif
hardening = ( p%h0 + p%h0_slopeLnRate * log(gamma_dot) ) &
* abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**p%a &
hardening = ( prm%h0 + prm%h0_slopeLnRate * log(gamma_dot) ) &
* abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**prm%a &
* sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation)
else
hardening = 0.0_pReal
@ -584,6 +498,7 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use material, only: &
plasticState, &
material_phase, &
phasememberAt, &
phase_plasticityInstance
@ -596,9 +511,9 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
ip, & !< integration point
el !< element
type(tParameters), pointer :: p
type(tParameters), pointer :: prm
real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults) :: &
plastic_isotropic_postResults
real(pReal), dimension(6) :: &
@ -613,11 +528,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)
prm => param(instance)
!--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress
if (p%dilatation) then
if (prm%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
@ -629,15 +544,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(p%outputID(o))
select case(prm%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) = &
p%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
prm%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
/ &!----------------------------------------------------------------------------------
(p%fTaylor * state(instance)%flowstress(of)) ) ** p%n
(prm%fTaylor * state(instance)%flowstress(of)) ) ** prm%n
c = c + 1_pInt
end select
enddo outputsLoop

49
src/plastic_kinematichardening.f90
View File

@ -1,6 +1,7 @@
!--------------------------------------------------------------------------------------------------
!> @author Philip Eisenlohr, Michigan State University
!> @author Zhuowen Zhao, Michigan State University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Introducing Voce-type kinematic hardening rule into crystal plasticity
!! formulation using a power law fitting
!--------------------------------------------------------------------------------------------------
@ -145,7 +146,8 @@ subroutine plastic_kinehardening_init(fileUnit)
phase_plasticityInstance, &
phase_Noutput, &
material_phase, &
plasticState, &
plasticState
use config, only: &
MATERIAL_partPhase
use lattice
use numerics,only: &
@ -155,9 +157,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 +180,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 +221,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)
@ -232,43 +232,53 @@ subroutine plastic_kinehardening_init(fileUnit)
allocate(param(instance)%nonSchmidCoeff(Nchunks_nonSchmid), source=0.0_pReal)
if(allocated(tempPerSlip)) deallocate(tempPerSlip)
allocate(tempPerSlip(Nchunks_SlipFamilies))
allocate(param(instance)%outputID(0))
endif
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 +629,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 +748,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 +806,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 +822,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 +874,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 +895,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)))) :: &

2
src/plastic_nonlocal.f90
View File

@ -291,8 +291,8 @@ use material, only: phase_plasticity, &
PLASTICITY_NONLOCAL_label, &
PLASTICITY_NONLOCAL_ID, &
plasticState, &
MATERIAL_partPhase ,&
material_phase
use config, only: MATERIAL_partPhase
use lattice
use numerics,only: &
numerics_integrator

3
src/plastic_phenopowerlaw.f90
View File

@ -157,7 +157,8 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
PLASTICITY_PHENOPOWERLAW_label, &
PLASTICITY_PHENOPOWERLAW_ID, &
material_phase, &
plasticState, &
plasticState
use config, only: &
MATERIAL_partPhase
use lattice
use numerics,only: &

1
src/porosity_none.f90
View File

@ -27,6 +27,7 @@ subroutine porosity_none_init()
use IO, only: &
IO_timeStamp
use material
use config
implicit none
integer(pInt) :: &

7
src/porosity_phasefield.f90
View File

@ -77,11 +77,10 @@ subroutine porosity_phasefield_init(fileUnit)
porosityState, &
porosityMapping, &
porosity, &
porosity_initialPhi, &
porosity_initialPhi
use config, only: &
material_partHomogenization, &
material_partPhase
use numerics,only: &
worldrank
implicit none
integer(pInt), intent(in) :: fileUnit
@ -94,11 +93,9 @@ subroutine porosity_phasefield_init(fileUnit)
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- porosity_'//POROSITY_phasefield_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(porosity_type == POROSITY_phasefield_ID),pInt)
if (maxNinstance == 0_pInt) return

5
src/source_damage_anisoBrittle.f90
View File

@ -91,9 +91,10 @@ subroutine source_damage_anisoBrittle_init(fileUnit)
phase_Noutput, &
SOURCE_damage_anisoBrittle_label, &
SOURCE_damage_anisoBrittle_ID, &
material_Nphase, &
material_phase, &
sourceState, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator

5
src/source_damage_anisoDuctile.f90
View File

@ -95,9 +95,10 @@ subroutine source_damage_anisoDuctile_init(fileUnit)
phase_Noutput, &
SOURCE_damage_anisoDuctile_label, &
SOURCE_damage_anisoDuctile_ID, &
material_Nphase, &
material_phase, &
sourceState, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator

5
src/source_damage_isoBrittle.f90
View File

@ -81,9 +81,10 @@ subroutine source_damage_isoBrittle_init(fileUnit)
phase_Noutput, &
SOURCE_damage_isoBrittle_label, &
SOURCE_damage_isoBrittle_ID, &
material_Nphase, &
material_phase, &
sourceState, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator

6
src/source_damage_isoDuctile.f90
View File

@ -81,10 +81,12 @@ subroutine source_damage_isoDuctile_init(fileUnit)
phase_Noutput, &
SOURCE_damage_isoDuctile_label, &
SOURCE_damage_isoDuctile_ID, &
material_Nphase, &
material_phase, &
sourceState, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator

5
src/source_thermal_dissipation.f90
View File

@ -67,9 +67,10 @@ subroutine source_thermal_dissipation_init(fileUnit)
phase_Noutput, &
SOURCE_thermal_dissipation_label, &
SOURCE_thermal_dissipation_ID, &
material_Nphase, &
material_phase, &
sourceState, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator

5
src/source_thermal_externalheat.f90
View File

@ -73,9 +73,10 @@ subroutine source_thermal_externalheat_init(fileUnit)
phase_Noutput, &
SOURCE_thermal_externalheat_label, &
SOURCE_thermal_externalheat_ID, &
material_Nphase, &
material_phase, &
sourceState, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator

5
src/source_vacancy_irradiation.f90
View File

@ -69,9 +69,10 @@ subroutine source_vacancy_irradiation_init(fileUnit)
phase_Noutput, &
SOURCE_vacancy_irradiation_label, &
SOURCE_vacancy_irradiation_ID, &
material_Nphase, &
material_phase, &
sourceState, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator

5
src/source_vacancy_phenoplasticity.f90
View File

@ -67,9 +67,10 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit)
phase_Noutput, &
SOURCE_vacancy_phenoplasticity_label, &
SOURCE_vacancy_phenoplasticity_ID, &
material_Nphase, &
material_phase, &
sourceState, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator

5
src/source_vacancy_thermalfluc.f90
View File

@ -71,9 +71,10 @@ subroutine source_vacancy_thermalfluc_init(fileUnit)
phase_Noutput, &
SOURCE_vacancy_thermalfluc_label, &
SOURCE_vacancy_thermalfluc_ID, &
material_Nphase, &
material_phase, &
sourceState, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator

46
src/spectral_mech_Basic.f90
View File

@ -2,7 +2,7 @@
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Basic scheme PETSc solver
!> @brief Basic scheme solver
!--------------------------------------------------------------------------------------------------
module spectral_mech_basic
#include <petsc/finclude/petscsnes.h>
@ -22,7 +22,7 @@ module spectral_mech_basic
private
character (len=*), parameter, public :: &
DAMASK_spectral_SolverBasicPETSC_label = 'basic'
DAMASK_spectral_SolverBasic_label = 'basic'
!--------------------------------------------------------------------------------------------------
! derived types
@ -65,9 +65,9 @@ module spectral_mech_basic
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
public :: &
basicPETSc_init, &
basicPETSc_solution, &
BasicPETSc_forward
basic_init, &
basic_solution, &
basic_forward
external :: &
PETScErrorF ! is called in the CHKERRQ macro
@ -76,7 +76,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine basicPETSc_init
subroutine basic_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
@ -124,9 +124,9 @@ subroutine basicPETSc_init
external :: &
SNESSetOptionsPrefix, &
SNESSetConvergenceTest, &
DMDASNESsetFunctionLocal
DMDASNESSetFunctionLocal
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:31–45, 2015'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -158,9 +158,9 @@ subroutine basicPETSc_init
call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
call DMsetUp(da,ierr); CHKERRQ(ierr)
call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,Basic_formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESsetConvergenceTest(snes,BasicPETSC_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
call SNESsetConvergenceTest(snes,Basic_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
CHKERRQ(ierr)
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
@ -212,12 +212,12 @@ subroutine basicPETSc_init
call Utilities_updateGamma(C_minMaxAvg,.true.)
end subroutine basicPETSc_init
end subroutine basic_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Basic PETSC scheme with internal iterations
!> @brief solution for the Basic scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
type(tSolutionState) function basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
use IO, only: &
IO_error
use numerics, only: &
@ -275,19 +275,19 @@ type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old,
! check convergence
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
BasicPETSc_solution%converged = reason > 0
basicPETSC_solution%iterationsNeeded = totalIter
basicPETSc_solution%termIll = terminallyIll
basic_solution%converged = reason > 0
basic_solution%iterationsNeeded = totalIter
basic_solution%termIll = terminallyIll
terminallyIll = .false.
if (reason == -4) call IO_error(893_pInt) ! MPI error
end function BasicPETSc_solution
end function basic_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the basic residual vector
!--------------------------------------------------------------------------------------------------
subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
subroutine Basic_formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: &
itmax, &
itmin
@ -370,13 +370,13 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
! constructing residual
f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
end subroutine BasicPETSc_formResidual
end subroutine Basic_formResidual
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
subroutine Basic_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
@ -425,14 +425,14 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end subroutine BasicPETSc_converged
end subroutine Basic_converged
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
subroutine Basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
@ -538,6 +538,6 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine BasicPETSc_forward
end subroutine Basic_forward
end module spectral_mech_basic

5
src/thermal_adiabatic.f90
View File

@ -64,6 +64,8 @@ subroutine thermal_adiabatic_init(fileUnit)
IO_error, &
IO_timeStamp, &
IO_EOF
use config, only: &
material_partHomogenization
use material, only: &
thermal_type, &
thermal_typeInstance, &
@ -76,8 +78,7 @@ subroutine thermal_adiabatic_init(fileUnit)
thermalMapping, &
thermal_initialT, &
temperature, &
temperatureRate, &
material_partHomogenization
temperatureRate
implicit none
integer(pInt), intent(in) :: fileUnit

3
src/thermal_conduction.f90
View File

@ -77,7 +77,8 @@ subroutine thermal_conduction_init(fileUnit)
thermalMapping, &
thermal_initialT, &
temperature, &
temperatureRate, &
temperatureRate
use config, only: &
material_partHomogenization
implicit none

1
src/thermal_isothermal.f90
View File

@ -27,6 +27,7 @@ subroutine thermal_isothermal_init()
use IO, only: &
IO_timeStamp
use material
use config
implicit none
integer(pInt) :: &

7
src/vacancyflux_cahnhilliard.f90
View File

@ -91,9 +91,10 @@ subroutine vacancyflux_cahnhilliard_init(fileUnit)
vacancyfluxMapping, &
vacancyConc, &
vacancyConcRate, &
vacancyflux_initialCv, &
material_partHomogenization, &
material_partPhase
vacancyflux_initialCv
use config, only: &
material_partPhase, &
material_partHomogenization
implicit none
integer(pInt), intent(in) :: fileUnit

3
src/vacancyflux_isochempot.f90
View File

@ -74,7 +74,8 @@ subroutine vacancyflux_isochempot_init(fileUnit)
vacancyfluxMapping, &
vacancyConc, &
vacancyConcRate, &
vacancyflux_initialCv, &
vacancyflux_initialCv
use config, only: &
material_partHomogenization
implicit none

1
src/vacancyflux_isoconc.f90
View File

@ -27,6 +27,7 @@ subroutine vacancyflux_isoconc_init()
use IO, only: &
IO_timeStamp
use material
use config
implicit none
integer(pInt) :: &