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

Merge remote-tracking branch 'origin/development' into homogenization-respect-fraction

This commit is contained in:
Martin Diehl 2022-02-14 09:19:49 +01:00
parent 28eda64893 1ecbeb6924
commit ebcf0fa625
92 changed files with 717 additions and 312 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
.gitlab-ci.yml
View File

@ -36,21 +36,21 @@ variables:
# Names of module files to load
# ===============================================================================================
# ++++++++++++ Compiler +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
COMPILER_GNU: "Compiler/GNU/10"
COMPILER_GNU: "Compiler/GNU/10"
COMPILER_INTELLLVM: "Compiler/oneAPI/2022.0.1 Libraries/IMKL/2022.0.1"
COMPILER_INTEL: "Compiler/Intel/2022.0.1 Libraries/IMKL/2022.0.1"
# ++++++++++++ MPI ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.1"
MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0"
MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.2"
MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0"
MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
# ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
PETSC_GNU: "Libraries/PETSc/3.16.4/GNU-10-OpenMPI-4.1.2"
PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.3/Intel-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.4/Intel-2022.0.1-IntelMPI-2021.5.0"
# ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MSC: "FEM/MSC/2021.3.1"
IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
HDF5Marc: "HDF5/1.12.1/Intel-19.1.2"
MSC: "FEM/MSC/2021.3.1"
IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
HDF5Marc: "HDF5/1.12.1/Intel-19.1.2"
###################################################################################################

2
CMakeLists.txt
View File

@ -42,7 +42,7 @@ string(TOUPPER "${CMAKE_BUILD_TYPE}" CMAKE_BUILD_TYPE)
if(CMAKE_BUILD_TYPE STREQUAL "DEBUG" OR CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
set(DEBUG_FLAGS "${DEBUG_FLAGS} -DDEBUG")
set(PARALLEL "OFF")
set(OPTI "OFF")
set(OPTI "DEBUG")
elseif(CMAKE_BUILD_TYPE STREQUAL "RELEASE")
set(PARALLEL "ON")
set(OPTI "DEFENSIVE")

2
PRIVATE

@ -1 +1 @@
Subproject commit 5774122bf48d637704bb4afb10b87c34a4dbcaba
Subproject commit 0d639a9ba41db279b0d2825c8e8eddf0ccd91326

17
cmake/Compiler-GNU.cmake
View File

@ -9,26 +9,24 @@ if (OPENMP)
set (OPENMP_FLAGS "-fopenmp")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
if (OPTIMIZATION STREQUAL "DEBUG")
set (OPTIMIZATION_FLAGS "-Og")
elseif (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2 -mtune=generic -flto")
set (OPTIMIZATION_FLAGS "-O2 -mtune=native")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-O3 -march=native -mtune=native -ffast-math -funroll-loops -ftree-vectorize -flto")
set (OPTIMIZATION_FLAGS "-O3 -march=native -funroll-loops -ftree-vectorize -flto")
endif ()
set (STANDARD_CHECK "-std=f2018 -pedantic-errors" )
set (LINKER_FLAGS "${LINKER_FLAGS} -Wl")
# options parsed directly to the linker
set (LINKER_FLAGS "${LINKER_FLAGS},-undefined,dynamic_lookup" )
# ensure to link against dynamic libraries
#------------------------------------------------------------------------------------------------
# Fine tuning compilation options
set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp")
# preprocessor
set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIC -fPIE")
set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIE")
# position independent code
set (COMPILE_FLAGS "${COMPILE_FLAGS} -ffree-line-length-132")
@ -123,6 +121,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -ffpe-trap=invalid,zero,overflow")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -g")
# Generate symbolic debugging information in the object file
set (DEBUG_FLAGS "${DEBUG_FLAGS} -Og")
# Optimize debugging experience
set (DEBUG_FLAGS "${DEBUG_FLAGS} -fbacktrace")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -fdump-core")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -fcheck=all")

7
cmake/Compiler-Intel.cmake
View File

@ -9,12 +9,12 @@ if (OPENMP)
set (OPENMP_FLAGS "-qopenmp -parallel")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG")
set (OPTIMIZATION_FLAGS "-O0 -no-ip")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost")
set (OPTIMIZATION_FLAGS "-ipo -O3 -fp-model fast=2 -xHost")
# -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost"
endif ()
@ -110,6 +110,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0")
# generate debug information for parameters
# Disabled due to ICE when compiling phase_damage.f90 (not understandable, there is no parameter in there)
set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
# generate complete debugging information
# Additional options
# -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
# -check: Checks at runtime, where

5
cmake/Compiler-IntelLLVM.cmake
View File

@ -9,7 +9,7 @@ if (OPENMP)
set (OPENMP_FLAGS "-qopenmp")
endif ()
if (OPTIMIZATION STREQUAL "OFF")
if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG")
set (OPTIMIZATION_FLAGS "-O0")
elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2")
@ -109,6 +109,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0")
set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all")
# generate debug information for parameters
set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
# generate complete debugging information
# Additional options
# -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
# -check: Checks at runtime, where

5
examples/config/phase/AISI304.yaml
View File

@ -1,6 +1,7 @@
references:
- H.M. Ledbetter
physica status solidi (a) 85(1):89-96, 1984
- H.M. Ledbetter,
physica status solidi (a) 85(1):89-96, 1984,
https://doi.org/10.1002/pssa.2210850111
lattice: cF
rho: 7937.0

1
examples/config/phase/Ag.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Silver
lattice: cF
rho: 10490.0

1
examples/config/phase/Al.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Aluminium
lattice: cF
rho: 2700.0

1
examples/config/phase/Au.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Gold
lattice: cF
rho: 19300.0

1
examples/config/phase/Cu.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Copper
lattice: cF
rho: 8960.0

1
examples/config/phase/Fe.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Iron
lattice: cI
rho: 7874.0

1
examples/config/phase/Mg.yaml
View File

@ -2,6 +2,7 @@ references:
- D. Tromans,
International Journal of Recent Research and Applied Studies 6(4):462-483, 2011,
https://www.arpapress.com/Volumes/Vol6Issue4/IJRRAS_6_4_14.pdf
lattice: hP
c/a: 1.62350
rho: 1740.0

1
examples/config/phase/Ni.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Nickel
lattice: cF
rho: 8908.0

9
examples/config/phase/Sn-beta.yaml Normal file
View File

@ -0,0 +1,9 @@
references:
- J.A. Rayne and B.S. Chandrasekhar,
Physical Review 120(5):1658-1663, 1960,
https://doi.org/10.1103/PhysRev.120.1658
- https://en.wikipedia.org/wiki/Tin
lattice: tI
c/a: 0.5458 # T=300K (c=31.83nm, a=5.832nm)
rho: 7265.0

1
examples/config/phase/Ti.yaml
View File

@ -1,6 +1,7 @@
references:
- https://www.totalmateria.com/page.aspx?ID=CheckArticle&site=ktn&NM=221
- https://en.wikipedia.org/wiki/Titanium
lattice: hP
c/a: 1.587
rho: 4506.0

1
examples/config/phase/W.yaml
View File

@ -1,4 +1,5 @@
references:
- https://en.wikipedia.org/wiki/Tungsten
lattice: cI
rho: 19300.0

2
examples/config/phase/bct.yaml Normal file
View File

@ -0,0 +1,2 @@
lattice: tI
c/a: 0.55

6
examples/config/phase/damage/anisobrittle_cubic.yaml
View File

@ -1,11 +1,13 @@
type: anisobrittle
output: [f_phi]
N_cl: [3]
g_crit: [0.5e+7]
s_crit: [0.006666]
dot_o: 1.e-3
q: 20
output: [f_phi]
D_11: 1.0
mu: 0.001

2
examples/config/phase/damage/isobrittle_generic.yaml
View File

@ -1,7 +1,7 @@
type: isobrittle
W_crit: 1400000.0
output: [f_phi]
W_crit: 1400000.0
D_11: 1.0
mu: 0.001

2
examples/config/phase/hcp.yaml Normal file
View File

@ -0,0 +1,2 @@
lattice: hP
c/a: 1.6333

16
examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml
View File

@ -1,7 +1,13 @@
type: thermalexpansion
references:
- R.H. Bogaard et al.
Thermochimica Acta 218:373-393, 1993
https://doi.org/10.1016/0040-6031(93)80437-F
A_11: 15.0e-6
T_ref: 300.0
- R.H. Bogaard et al.,
Thermochimica Acta 218:373-393, 1993,
https://doi.org/10.1016/0040-6031(93)80437-F,
fit to Fig. 6 (T_min=100K, T_max=1400K)
A_11: 2.068e-08
A_11,T: 1.579e-09
A_11,T^2: 3.449e-13
T_ref: 293.15

6
examples/config/phase/mechanical/eigen/thermalexpansion_Al.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 23.1e-6
T_ref: 293.15

6
examples/config/phase/mechanical/eigen/thermalexpansion_Au.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 14.e-6
T_ref: 293.15

5
examples/config/phase/mechanical/eigen/thermalexpansion_C35E.yaml
View File

@ -1,8 +1,11 @@
type: thermalexpansion
references:
- https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg,
fitted from image description (Scilab code)
fit to image description (Scilab code)
A_11: 12.70371e-6
A_11,T: 7.54e-9
A_11,T^2: -1.0e-11
T_ref: 273.0

6
examples/config/phase/mechanical/eigen/thermalexpansion_Cu.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 17.e-6
T_ref: 293.15

6
examples/config/phase/mechanical/eigen/thermalexpansion_Fe.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 11.8e-6
T_ref: 293.15

17
examples/config/phase/mechanical/eigen/thermalexpansion_Sn-beta.yaml Normal file
View File

@ -0,0 +1,17 @@
type: thermalexpansion
references:
- V.T. Deshpande and D.B. Sirdeshmukh,
Acta Crystallographica 15:294-295, 1962,
https://doi.org/10.1107/S0365110X62000742,
fit to Tab. 2 (T_min=30ºC, T_max=210ºC)
A_11: 1.639e-05
A_11,T: 1.799e-08
A_11,T^2: 1.734e-10
A_33: 3.263e-05
A_33,T: 1.387e-08
A_33,T^2: 5.794e-10
T_ref: 293.15

6
examples/config/phase/mechanical/eigen/thermalexpansion_W.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion
references:
- https://en.wikipedia.org/wiki/Thermal_expansion
- https://en.wikipedia.org/wiki/Thermal_expansion,
293.15K
A_11: 4.5e-6
T_ref: 293.15

7
examples/config/phase/mechanical/eigen/thermalexpansion_X20Cr13.yaml
View File

@ -1,7 +1,10 @@
type: thermalexpansion
references:
- https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg
fitted from image description (Scilab code)
- https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg,
fit to image description (Scilab code)
A_11: 11.365e-6
A_11,T: 5.0e-9
T_ref: 273.0

6
examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml
View File

@ -1,8 +1,10 @@
type: Hooke
references:
- H.M. Ledbetter
physica status solidi (a) 85(1):89-96, 1984
- H.M. Ledbetter,
physica status solidi (a) 85(1):89-96, 1984,
https://doi.org/10.1002/pssa.2210850111
C_11: 204.6e+9
C_12: 137.7e+9
C_44: 126.2e+9

5
examples/config/phase/mechanical/elastic/Hooke_Ag.yaml
View File

@ -1,4 +1,5 @@
type: Hooke
references:
- J.R. Neighbours and G.A. Alers,
Physical Review 111:707-712, 1958,
@ -7,8 +8,6 @@ references:
Journal of Applied Physics 37:3567-3572, 1966,
https://doi.org/10.1063/1.1708903
T_ref: 300
C_11: 122.9e+9
C_11,T: -313.5e+5
C_11,T^2: -107.3e+2
@ -20,3 +19,5 @@ C_12,T^2: -681.6e+1
C_44: 42.63e+9
C_44,T: -180.5e+5
C_44,T^2: -353.8e+1
T_ref: 300

27
examples/config/phase/mechanical/elastic/Hooke_Al.yaml
View File

@ -1,22 +1,25 @@
type: Hooke
references:
- G.N. Kamm and G.A. Alers,
Journal of Applied Physics 35:327-330, 1964,
https://doi.org/10.1063/1.1713309
https://doi.org/10.1063/1.1713309,
fit to Tab. I (T_min=100K, T_max=300K)
- D. Gerlich and E.S. Fisher,
Journal of Physics and Chemistry of Solids 30:1197-1205, 1969
https://doi.org/10.1016/0022-3697(69)90377-1
https://doi.org/10.1016/0022-3697(69)90377-1,
fit to Tab. 2 (T_min=293K, T_max=900K)
T_ref: 300
C_11: 106.9e+9
C_11,T: -3.685e+7
C_11,T^2: -1.016e+4
C_11: 106.1e+9
C_11,T: -359.3e+5
C_11,T^2: -152.7e+2
C_12: 60.55e+9
C_12,T: -8.307e+6
C_12,T^2: -4.353e+3
C_12: 57.83e+9
C_12,T: -781.6e+4
C_12,T^2: -551.3e+1
C_44: 28.37e+9
C_44,T: -1.418e+7
C_44,T^2: -3.253e+3
C_44: 24.31e+9
C_44,T: -142.9e+5
C_44,T^2: -404.6e+1
T_ref: 293.15

2
examples/config/phase/mechanical/elastic/Hooke_Au.yaml
View File

@ -1,9 +1,11 @@
type: Hooke
references:
- J.P. Hirth and J. Lothe,
Theory of Dislocations, 1982,
John Wiley & Sons,
page 837
C_11: 186.e+9
C_12: 157.e+9
C_44: 42.e+9

24
examples/config/phase/mechanical/elastic/Hooke_Cu.yaml
View File

@ -1,7 +1,21 @@
type: Hooke
references:
- https://www.mit.edu/~6.777/matprops/copper.htm,
fixed typo
C_11: 168.3e+9
C_12: 122.1e+9
C_44: 75.7e+9
- W.C. Overton, Jr. and J. Gaffney,
Physical Review 98(4):969-977, 1955,
https://doi.org/10.1103/PhysRev.98.969,
fit to Tab. I (T_min=100K, T_max=300K)
C_11: 168.6e+9
C_11,T: -3.779e+7
C_11,T^2: -2.536e+4
C_12: 121.5e+9
C_12,T: -1.632e+7
C_12,T^2: -1.116e+4
C_44: 75.59e+9
C_44,T: -2.912e+7
C_44,T^2: -1.669e+4
T_ref: 293.15

22
examples/config/phase/mechanical/elastic/Hooke_Fe.yaml
View File

@ -1,19 +1,21 @@
type: Hooke
references:
- D.J. Dever,
Journal of Applied Physics 43(8):3293-3301, 1972,
https://doi.org/10.1063/1.1661710
fit to Tab. II (T_min=25ºC, T_max=880ºC)
T_ref: 300
C_11: 232.1e+9
C_11,T: -4.678e+7
C_11,T^2: -5.762e+4
C_11: 231.7e+9
C_11,T: -47.6e+6
C_11,T^2: -54.4e+3
C_12: 135.9e+9
C_12,T: -1.695e+7
C_12,T^2: 1.555e+3
C_12: 135.8e+9
C_12,T: -12.9e+6
C_12,T^2: -7.3e+3
C_44: 117.0e+9
C_44,T: -2.047e+7
C_44,T^2: -2.814e+2
C_44: 116.8e+9
C_44,T: -19.4e+6
C_44,T^2: -2.5e+3
T_ref: 293.15

35
examples/config/phase/mechanical/elastic/Hooke_Mg.yaml
View File

@ -1,10 +1,29 @@
type: Hooke
references:
- D. Tromans,
International Journal of Recent Research and Applied Studies 6(4):462-483, 2011,
https://www.arpapress.com/Volumes/Vol6Issue4/IJRRAS_6_4_14.pdf
C_11: 59.3e+9
C_33: 61.5e+9
C_44: 16.4e+9
C_12: 25.7e+9
C_13: 21.4e+9
- L.J. Slutsky and C.W. Garland,
Physical Review 107(4):972-976, 1957,
https://doi.org/10.1103/PhysRev.107.972,
fit to Tab. I (T_min=100K, T_max=300K)
C_11: 59.50e+9
C_11,T: -1.930e+7
C_11,T^2: -1.215e+4
C_33: 61.72e+9
C_33,T: -2.175e+7
C_33,T^2: -5.624e+3
C_44: 16.46e+9
C_44,T: -1.006e+7
C_44,T^2: -7.692e+3
C_12: 25.62e+9
C_12,T: -2.216e+6
C_12,T^2: -4.138e+3
C_13: 21.46e+9
C_13,T: -1.921e+6
C_13,T^2: -4.283e+3
T_ref: 293.15

26
examples/config/phase/mechanical/elastic/Hooke_Ni.yaml
View File

@ -1,9 +1,21 @@
type: Hooke
references:
- J.P. Hirth and J. Lothe,
Theory of Dislocations, 1982,
John Wiley & Sons,
page 837
C_11: 246.5e+9
C_12: 147.3e+9
C_44: 124.7e+9
- G.A. Alers,
Journal of Physics and Chemistry of Solids 13(1-2):40-55, 1960,
https://doi.org/10.1016/0022-3697(60)90125-6,
fit to Tab. 2 (T_min=100K, T_max=700K)
C_11: 251.0e+9
C_11,T: -4.998e+7
C_11,T^2: -2.952e+4
C_12: 150.0e+9
C_12,T: -4.269e+6
C_12,T^2: -6.429e+3
C_44: 123.7e+9
C_44,T: -3.618e+7
C_44,T^2: -7.024e+3
T_ref: 293.15

2
examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml
View File

@ -1,8 +1,10 @@
type: Hooke
references:
- S.A. Kim and W.L. Johnson,
Materials Science & Engineering A 452-453:633-639, 2007,
https://doi.org/10.1016/j.msea.2006.11.147
C_11: 268.1e+9
C_12: 111.2e+9
C_44: 79.06e+9

32
examples/config/phase/mechanical/elastic/Hooke_Sn-beta.yaml Normal file
View File

@ -0,0 +1,32 @@
type: Hooke
references:
- J.A. Rayne and B.S. Chandrasekhar,
Physical Review 120(5):1658-1663, 1960,
https://doi.org/10.1103/PhysRev.120.1658,
fit to Fig. 2 (T_min=100K, T_max=300K) and Tab. IV (C_13, T_min=77K, T_max=300K)
C_11: 72.90e+9
C_11,T: -4.399e+7
C_11,T^2: -2.645e+4
C_12: 59.27e+9
C_12,T: 1.058e+7
C_12,T^2: 1.002e+4
C_13: 35.80e+9
C_13,T: -2.870e+6
C_33: 88.78e+9
C_33,T: -5.250e+7
C_33,T^2: 3.546e+3
C_44: 22.26e+9
C_44,T: -1.982e+7
C_44,T^2: -8.711e+3
C_66: 24.18e+9
C_66,T: -1.806e+7
C_66,T^2: -4.112e+3
T_ref: 293.15

2
examples/config/phase/mechanical/elastic/Hooke_TWIP-steel.yaml
View File

@ -1,4 +1,5 @@
type: Hooke
references:
- D. Music et al.,
Applied Physics Letters 99(19):191904, 2007,
@ -6,6 +7,7 @@ references:
- S.L. Wong et al.,
Acta Materialia 118:140-151, 2016,
https://doi.org/10.1016/j.actamat.2016.07.032
C_11: 175.0e+9
C_12: 115.0e+9
C_44: 135.0e+9

39
examples/config/phase/mechanical/elastic/Hooke_Ti.yaml
View File

@ -1,10 +1,33 @@
type: Hooke
references:
- L. Wang et al.,
Acta Materialia 132:598-610, 2017,
https://doi.org/10.1016/j.actamat.2017.05.015
C_11: 162.4e+9
C_33: 181.6e+9
C_44: 47.2e+9
C_12: 92.e+9
C_13: 69.e+9
- E.S. Fisher and C.J. Renken,
Physical Review 135(2A):A482-A494, 1964,
https://doi.org/10.1103/PhysRev.135.A482,
fit to Tab. IV (T_min=150K, T_max=250K)
- H. Ogi et al.,
Acta Materialia 52(7):2075-2080, 2004,
https://doi.org/10.1016/j.actamat.2004.01.002,
fit to Fig. 3 (T_min=300K, T_max=900K)
C_11: 162.6e+9
C_11,T: -6.150e+7
C_11,T^2: -5.557e+2
C_33: 183.3e+9
C_33,T: -1.655e+07
C_33,T^2: -1.022e+04
C_44: 45.80e+9
C_44,T: -2.936e+07
C_44,T^2: 7.120e+02
C_12: 89.97e+9
C_12,T: 2.776e+6
C_12,T^2: -2.389e+4
C_13: 69.53e+9
C_13,T: 1.057e+7
C_13,T^2: -2.966e+3
T_ref: 293.15

24
examples/config/phase/mechanical/elastic/Hooke_W.yaml
View File

@ -1,8 +1,20 @@
type: Hooke
references:
- D. Cereceda et al.,
International Journal of Plasticity 78:242-265, 2016,
https://doi.org/10.1016/j.ijplas.2015.09.002
C_11: 523.e+9
C_12: 202.e+9
C_44: 161.e+9
- F.H. Featherston and J.R. Nieghbours,
Physical Review 130(4):1324-1333,
https://doi.org/10.1103/PhysRev.130.1324,
fit to Tab. III (T_min=100K, T_max=300K)
C_11: 523.6e+9
C_11,T: -7.607e+7
C_11,T^2: -1.551e+5
C_12: 205.1e+9
C_12,T: -2.843e+6
C_44: 160.8e+9
C_44,T: -1.057e+7
C_44,T^2: 9.933e+3
T_ref: 293.15

2
examples/config/phase/mechanical/elastic/Hooke_vanishing-Poisson-ratio.yaml
View File

@ -1,8 +1,10 @@
type: Hooke
references:
- T. Maiti and P. Eisenlohr,
Scripta Materialia 145:37-40, 2018,
https://doi.org/10.1016/j.scriptamat.2017.09.047
C_11: 1.e+8
C_12: 1.e+6
C_44: 4.95e+7

4
examples/config/phase/mechanical/plastic/dislotungsten_W.yaml
View File

@ -1,4 +1,5 @@
type: dislotungsten
references:
- D. Cereceda et al.,
International Journal of Plasticity 78:242-265, 2016,
@ -6,8 +7,11 @@ references:
- R. Gröger et al.,
Acta Materialia 56(19):5412-5425, 2008,
https://doi.org/10.1016/j.actamat.2008.07.037
output: [Lambda_sl]
N_sl: [12]
b_sl: [2.72e-10]
rho_mob_0: [1.0e+9] # estimated from section 3.2
rho_dip_0: [1.0] # not given

4
examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml
View File

@ -1,4 +1,5 @@
type: dislotwin
references:
- K. Sedighiani et al.,
International Journal of Plasticity 134:102779, 2020,
@ -6,8 +7,11 @@ references:
- K. Sedighiani et al.,
Mechanics of Materials, 164:104117, 2022,
https://doi.org/10.1016/j.mechmat.2021.104117
output: [rho_dip, rho_mob]
N_sl: [12, 12]
b_sl: [2.49e-10, 2.49e-10]
rho_mob_0: [2.81e12, 2.8e+12]
rho_dip_0: [1.0, 1.0] # not given

3
examples/config/phase/mechanical/plastic/isotropic_free-surface.yaml
View File

@ -1,9 +1,12 @@
type: isotropic
references:
- T. Maiti and P. Eisenlohr,
Scripta Materialia 145:37-40, 2018,
https://doi.org/10.1016/j.scriptamat.2017.09.047
output: [xi]
dot_gamma_0: 0.001
n: 20.
xi_0: 0.3e+6

11
examples/config/phase/mechanical/plastic/nonlocal_Al.yaml
View File

@ -1,10 +1,13 @@
type: nonlocal
references:
C. Kords,
On the role of dislocation transport in the constitutive description of crystal plasticity,
RWTH Aachen 2013,
http://publications.rwth-aachen.de/record/229993/files/4862.pdf
- C. Kords,
On the role of dislocation transport in the constitutive description of crystal plasticity,
RWTH Aachen 2013,
http://publications.rwth-aachen.de/record/229993/files/4862.pdf
output: [rho_u_ed_pos, rho_b_ed_pos, rho_u_ed_neg, rho_b_ed_neg, rho_u_sc_pos, rho_b_sc_pos, rho_u_sc_neg, rho_b_sc_neg, rho_d_ed, rho_d_sc]
N_sl: [12]
b_sl: [2.86e-10]

11
examples/config/phase/mechanical/plastic/nonlocal_Ni.yaml
View File

@ -1,10 +1,13 @@
type: nonlocal
references:
C. Kords,
On the role of dislocation transport in the constitutive description of crystal plasticity,
RWTH Aachen 2013,
http://publications.rwth-aachen.de/record/229993/files/4862.pdf
- C. Kords,
On the role of dislocation transport in the constitutive description of crystal plasticity,
RWTH Aachen 2013,
http://publications.rwth-aachen.de/record/229993/files/4862.pdf
output: [rho_u_ed_pos, rho_b_ed_pos, rho_u_ed_neg, rho_b_ed_neg, rho_u_sc_pos, rho_b_sc_pos, rho_u_sc_neg, rho_b_sc_neg, rho_d_ed, rho_d_sc]
N_sl: [12]
b_sl: [2.48e-10]

4
examples/config/phase/mechanical/plastic/phenopowerlaw_Al.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- W.F. Hosford et al.,
Acta Metallurgica 8(3):187-199, 1960,
@ -7,8 +8,11 @@ references:
- U.F. Kocks,
Metallurgical and Materials Transactions B 1:1121–1143, 1970,
https://doi.org/10.1007/BF02900224
output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 20
a_sl: 3.1
h_0_sl-sl: 1.7e+8

4
examples/config/phase/mechanical/plastic/phenopowerlaw_Au.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- D. Ma et al.,
Acta Materialia 103:796-808, 2016,
@ -9,8 +10,11 @@ references:
- U.F. Kocks,
Metallurgical and Materials Transactions B 1:1121–1143, 1970,
https://doi.org/10.1007/BF02900224
output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 83.3
a_sl: 1.0
h_0_sl-sl: 75.0e+6

4
examples/config/phase/mechanical/plastic/phenopowerlaw_Cu.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- T Takeuchi,
Transactions of the Japan Institute of Metals 16(10):629-640, 1975,
@ -7,8 +8,11 @@ references:
- U.F. Kocks,
Metallurgical and Materials Transactions B 1:1121–1143, 1970,
https://doi.org/10.1007/BF02900224
output: [xi_sl, gamma_sl]
N_sl: [12]
n_sl: 20
a_sl: 1.0
h_0_sl-sl: 2.4e+8

4
examples/config/phase/mechanical/plastic/phenopowerlaw_DP-steel-ferrite.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- C.C. Tasan et al.,
Acta Materialia 81:386-400, 2014,
@ -6,8 +7,11 @@ references:
- U.F. Kocks,
Metallurgical and Materials Transactions B 1:1121–1143, 1970,
https://doi.org/10.1007/BF02900224
output: [xi_sl, gamma_sl]
N_sl: [12, 12]
n_sl: 20
a_sl: 2.25
h_0_sl-sl: 1.0e+9

1
examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- F. Wang et al.,
Acta Materialia 80:77-93, 2014,

33
examples/config/phase/mechanical/plastic/phenopowerlaw_Sn-beta.yaml Normal file
View File

@ -0,0 +1,33 @@
type: phenopowerlaw
references:
- A. Chakraborty and P. Eisenlohr,
Journal of Applied Physics 124:025302, 2018,
https://doi.org/10.1063/1.5029933
output: [xi_sl, gamma_sl]
N_sl: [2, 2, 2, 4, 2, 4, 2, 2, 4, 0, 0, 8]
n_sl: 6.0
a_sl: 2.0
h_0_sl-sl: 20.0e+6
xi_0_sl: [8.5e+6, 4.3e+6, 10.4e+6, 4.5e+6, 5.6e+6, 5.1e+6, 7.4e+6, 15.0e+6, 6.6e+6, 0.0, 0.0, 12.0e+6]
xi_inf_sl: [11.0e+6, 9.0e+6, 11.0e+6, 9.0e+6, 10.0e+6, 10.0e+6, 10.0e+6, 10.0e+6, 9.0e+6, 0.0, 0.0, 13.0e+6]
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, # 50
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, # 100
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, # 150
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # unused entries are indicated by -1.0
dot_gamma_0_sl: 2.6e-8

3
examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml
View File

@ -1,4 +1,5 @@
type: phenopowerlaw
references:
- C. Zambaldi et al.,
Journal of Materials Research 27(1):356-367, 2021,
@ -6,9 +7,11 @@ references:
- L. Wang et al.,
Acta Materialia 132:598-610, 2017,
https://doi.org/10.1016/j.actamat.2017.05.015
output: [gamma_sl]
N_sl: [3, 3, 0, 12] # basal, prism, -, 1. pyr<c+a>
n_sl: 20
a_sl: 2.0
dot_gamma_0_sl: 0.001

1
examples/config/phase/thermal/AISI304.yaml
View File

@ -5,5 +5,6 @@ references:
- R.H. Bogaard et al.
Thermochimica Acta 218:373-393, 1993
https://doi.org/10.1016/0040-6031(93)80437-F
C_p: 470.0
K_11: 14.34

13
examples/config/phase/thermal/Al.yaml
View File

@ -1,5 +1,14 @@
references:
- https://www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
- J.G. Hust and A.B. Lankford,
Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
US Department of Commerce, Boulder, Colorado, 1984,
fit to Tab. 3.4.1 (RRR=1000, T_min=200K, T_max=900K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
K_11: 2.380e+2
K_11,T: 2.068e-3
K_11,T^2: -7.765e-5
T_ref: 293.15
C_p: 910.0
K_11: 236.0

1
examples/config/phase/thermal/Au.yaml
View File

@ -1,4 +1,5 @@
references:
- https://de.wikipedia.org/wiki/Gold
C_p: 128.0
K_11: 320.0

12
examples/config/phase/thermal/Cu.yaml
View File

@ -1,4 +1,14 @@
references:
- J.G. Hust and A.B. Lankford,
Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
US Department of Commerce, Boulder, Colorado, 1984,
fit to Tab. 2.4.1 (RRR=1000, T_min=200K, T_max=1000K)
- https://www.mit.edu/~6.777/matprops/copper.htm
K_11: 4.039e+2
K_11,T: -8.119e-2
K_11,T^2: 1.454e-5
T_ref: 293.15
C_p: 385.0
K_11: 401.0

14
examples/config/phase/thermal/Fe.yaml Normal file
View File

@ -0,0 +1,14 @@
references:
- J.G. Hust and A.B. Lankford,
Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
US Department of Commerce, Boulder, Colorado, 1984,
fit to Tab. 4.4.1 (RRR=300, T_min=200K, T_max=1000K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
K_11: 8.055e+1
K_11,T: -1.051e-1
K_11,T^2: 5.464e-5
T_ref: 293.15
C_p: 450.0

14
examples/config/phase/thermal/Ni.yaml Normal file
View File

@ -0,0 +1,14 @@
references:
- Y.S. Touloukian et al.,
TPRC Data Series Volume 1. Thermal conductivity - metallic elements and alloys,
IFI/Plenum, 1970,
fit to Tab. 35R (T_min=150K, T_max=500K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
K_11: 9.132e+1
K_11,T: -1.525e-1
K_11,T^2: 3.053e-4
T_ref: 293.15
C_p: 440.0

18
examples/config/phase/thermal/Sn-beta.yaml Normal file
View File

@ -0,0 +1,18 @@
references:
- Y.S. Touloukian et al.,
TPRC Data Series Volume 1. Thermal conductivity - metallic elements and alloys,
IFI/Plenum, 1970,
fit to Tab. 61R (T_min=100K, T_max=400K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
K_11: 7.414e+1
K_11,T: -6.465e-2
K_11,T^2: 2.066e-4
K_33: 5.147e+1
K_33,T: -4.506e-2
K_33,T^2: 1.435e-4
T_ref: 293.15
C_p: 210.0

12
examples/config/phase/thermal/W.yaml
View File

@ -1,4 +1,14 @@
references:
- J.G. Hust and A.B. Lankford,
Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
US Department of Commerce, Boulder, Colorado, 1984,
fit to Tab. 5.4.1 (RRR=300, T_min=200K, T_max=1000K)
- https://www.mit.edu/~6.777/matprops/tungsten.htm
K_11: 1.758e+2
K_11,T: -1.605e-1
K_11,T^2: 1.160e-4
T_ref: 293.15
C_p: 132.51
K_11: 178.0

1
examples/config/phase/thermal/source/dissipation_generic.yaml
View File

@ -1,2 +1,3 @@
type: dissipation
kappa: .9

1
examples/config/phase/thermal/source/externalheat_ramp-and-hold.yaml
View File

@ -1,3 +1,4 @@
type: externalheat
f_T: [1, 1, 0, 0]
t_n: [0, 500, 500.001, 1000]

1
examples/config/phase/thermal/steel-0.5C.yaml
View File

@ -1,5 +1,6 @@
references:
- https://www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
C_p: 490.0
K_11: 54.0

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha5-556-g97f849c09
v3.0.0-alpha5-624-g26979da58

1
src/CPFEM.f90
View File

@ -90,7 +90,6 @@ subroutine CPFEM_initAll
call material_init(.false.)
call phase_init
call homogenization_init
call crystallite_init
call CPFEM_init
call config_deallocate

1
src/CPFEM2.f90
View File

@ -68,7 +68,6 @@ subroutine CPFEM_initAll
call material_init(restart=interface_restartInc>0)
call phase_init
call homogenization_init
call crystallite_init
call CPFEM_init
call config_deallocate

35
src/HDF5_utilities.f90
View File

@ -10,7 +10,8 @@ module HDF5_utilities
#include <petsc/finclude/petscsys.h>
use PETScSys
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
use MPI
use MPI_f08
use MPI, only: MPI_INFO_NULL_F90 => MPI_INFO_NULL
#endif
#endif
@ -162,9 +163,9 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel)
character, intent(in), optional :: mode
logical, intent(in), optional :: parallel
character :: m
integer(HID_T) :: plist_id
integer :: hdferr
character :: m
integer(HID_T) :: plist_id
integer :: hdferr
if (present(mode)) then
@ -178,9 +179,15 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel)
#ifdef PETSC
if (present(parallel)) then
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr)
else
call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr)
#else
if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
else
call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
#endif
end if
if(hdferr < 0) error stop 'HDF5 error'
#endif
@ -1860,7 +1867,7 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
globalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
integer, dimension(worldsize) :: &
integer(MPI_INTEGER_KIND), dimension(worldsize) :: &
readSize !< contribution of all processes
integer :: hdferr
integer(MPI_INTEGER_KIND) :: err_MPI
@ -1871,13 +1878,13 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
if(hdferr < 0) error stop 'HDF5 error'
!--------------------------------------------------------------------------------------------------
readSize = 0
readSize(worldrank+1) = int(localShape(ubound(localShape,1)))
readSize = 0_MPI_INTEGER_KIND
readSize(worldrank+1) = int(localShape(ubound(localShape,1)),MPI_INTEGER_KIND)
#ifdef PETSC
if (parallel) then
call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process
call MPI_Allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
end if
#endif
@ -1954,8 +1961,8 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
totalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id
integer, dimension(worldsize) :: writeSize !< contribution of all processes
integer(HID_T) :: dcpl
integer(MPI_INTEGER_KIND), dimension(worldsize) :: writeSize !< contribution of all processes
integer(HID_T) :: dcpl
integer :: hdferr
integer(MPI_INTEGER_KIND) :: err_MPI
integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T
@ -1974,11 +1981,11 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
!--------------------------------------------------------------------------------------------------
! determine the global data layout among all processes
writeSize = 0
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)))
writeSize = 0_MPI_INTEGER_KIND
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)),MPI_INTEGER_KIND)
#ifdef PETSC
if (parallel) then
call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process
call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
end if
#endif
@ -2009,7 +2016,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
if (hdferr < 0) error stop 'HDF5 error'
end if
end if
!--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape) and in file (global shape)
call H5Screate_simple_f(size(myShape), myShape, memspace_id, hdferr, myShape)

4
src/IO.f90
View File

@ -584,8 +584,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
character(len=pStringLen) :: formatString
select case (warning_ID)
case (42)
msg = 'parameter has no effect'
case (47)
msg = 'no valid parameter for FFTW, using FFTW_PATIENT'
case (207)
@ -594,8 +592,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
msg = 'crystallite responds elastically'
case (601)
msg = 'stiffness close to zero'
case (700)
msg = 'unknown crystal symmetry'
case (709)
msg = 'read only the first document'
case default

2
src/grid/DAMASK_grid.f90
View File

@ -32,7 +32,7 @@ program DAMASK_grid
implicit none
type :: tLoadCase
type(tRotation) :: rot !< rotation of BC
type(tRotation) :: rot !< rotation of BC
type(tBoundaryCondition) :: stress, & !< stress BC
deformation !< deformation BC (dot_F, F, or L)
real(pReal) :: t, & !< length of increment

2
src/grid/grid_mech_FEM.f90
View File

@ -339,7 +339,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(tRotation), intent(in) :: &
type(tRotation), intent(in) :: &
rotation_BC
PetscErrorCode :: err_PETSc
PetscScalar, pointer, dimension(:,:,:,:) :: &

8
src/grid/grid_mech_spectral_basic.f90
View File

@ -79,6 +79,12 @@ module grid_mechanical_spectral_basic
err_BC, & !< deviation from stress BC
err_div !< RMS of div of P
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
type(MPI_Status) :: status
#else
integer, dimension(MPI_STATUS_SIZE) :: status
#endif
integer :: &
totalIter = 0 !< total iteration in current increment
@ -244,7 +250,7 @@ subroutine grid_mechanical_spectral_basic_init
call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', &
MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI)
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_close(fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'

8
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -90,6 +90,12 @@ module grid_mechanical_spectral_polarisation
err_curl, & !< RMS of curl of F
err_div !< RMS of div of P
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
type(MPI_Status) :: status
#else
integer, dimension(MPI_STATUS_SIZE) :: status
#endif
integer :: &
totalIter = 0 !< total iteration in current increment
@ -270,7 +276,7 @@ subroutine grid_mechanical_spectral_polarisation_init
call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', &
MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI)
call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_File_close(fileUnit,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'

2
src/math.f90
View File

@ -898,7 +898,7 @@ end function math_33toVoigt6_stress
!--------------------------------------------------------------------------------------------------
!> @brief Convert 3x3 stress tensor into 6 Voigt vector.
!> @brief Convert 3x3 strain tensor into 6 Voigt vector.
!--------------------------------------------------------------------------------------------------
pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde)

4
src/mesh/DAMASK_mesh.f90
View File

@ -307,9 +307,11 @@ program DAMASK_mesh
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
end if
if (.not. cutBack .and. worldrank == 0) &
if (.not. cutBack .and. worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
flush(statUnit)
endif
end do subStepLooping
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc

20
src/parallelization.f90
View File

@ -52,13 +52,13 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine parallelization_init
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize, version, subversion, devNull
character(len=4) :: rank_str
character(len=MPI_MAX_LIBRARY_VERSION_STRING) :: MPI_library_version
!$ integer :: got_env, threadLevel
!$ integer(pI32) :: OMP_NUM_THREADS
!$ character(len=6) NumThreadsString
PetscErrorCode :: err_PETSc
#ifdef _OPENMP
! If openMP is enabled, check if the MPI libary supports it and initialize accordingly.
@ -86,12 +86,22 @@ subroutine parallelization_init
if (err_MPI /= 0_MPI_INTEGER_KIND) &
error stop 'Could not determine worldrank'
if (worldrank == 0) print'(/,1x,a)', '<<<+- parallelization init -+>>>'
if (worldrank == 0) then
print'(/,1x,a)', '<<<+- parallelization init -+>>>'
call MPI_Get_library_version(MPI_library_version,devNull,err_MPI)
print'(/,1x,a)', trim(MPI_library_version)
call MPI_Get_version(version,subversion,err_MPI)
print'(1x,a,i0,a,i0)', 'MPI standard: ',version,'.',subversion
#ifdef _OPENMP
print'(1x,a,i0)', 'OpenMP version: ',openmp_version
#endif
end if
call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) &
error stop 'Could not determine worldsize'
if (worldrank == 0) print'(/,1x,a,i3)', 'MPI processes: ',worldsize
if (worldrank == 0) print'(/,1x,a,i0)', 'MPI processes: ',worldsize
call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) &
@ -128,7 +138,7 @@ subroutine parallelization_init
!$ OMP_NUM_THREADS = 4_pI32
!$ endif
!$ endif
!$ print'(1x,a,1x,i2)', 'OMP_NUM_THREADS:',OMP_NUM_THREADS
!$ print'(1x,a,i0)', 'OMP_NUM_THREADS: ',OMP_NUM_THREADS
!$ call omp_set_num_threads(OMP_NUM_THREADS)
end subroutine parallelization_init

26
src/phase.f90
View File

@ -278,6 +278,7 @@ module phase
en
end subroutine plastic_dependentState
module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
integer, intent(in) :: ph, en
real(pReal), intent(in), dimension(3,3) :: &
@ -321,7 +322,6 @@ module phase
phase_restore, &
plastic_nonlocal_updateCompatibility, &
converged, &
crystallite_init, &
phase_mechanical_constitutive, &
phase_thermal_constitutive, &
phase_damage_constitutive, &
@ -399,6 +399,8 @@ subroutine phase_init
call damage_init
call thermal_init(phases)
call crystallite_init()
end subroutine phase_init
@ -434,6 +436,8 @@ subroutine phase_allocateState(state, &
allocate(state%dotState (sizeDotState,NEntries), source=0.0_pReal)
allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal)
state%deltaState2 => state%state(state%offsetDeltaState+1: &
state%offsetDeltaState+state%sizeDeltaState,:)
end subroutine phase_allocateState
@ -499,22 +503,12 @@ subroutine crystallite_init()
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el, & !< counter in element loop
cMax, & !< maximum number of integration point components
iMax, & !< maximum number of integration points
eMax, & !< maximum number of elements
en, ph
class(tNode), pointer :: &
num_crystallite, &
phases
print'(/,1x,a)', '<<<+- crystallite init -+>>>'
cMax = homogenization_maxNconstituents
iMax = discretization_nIPs
eMax = discretization_Nelems
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal)
@ -548,15 +542,9 @@ subroutine crystallite_init()
phases => config_material%get('phase')
print'(/,a42,1x,i10)', ' # of elements: ', eMax
print'( a42,1x,i10)', ' # of integration points/element: ', iMax
print'( a42,1x,i10)', 'max # of constituents/integration point: ', cMax
flush(IO_STDOUT)
!$OMP PARALLEL DO PRIVATE(ce,ph,en)
do el = 1, eMax
do ip = 1, iMax
do el = 1, discretization_Nelems
do ip = 1, discretization_nIPs
ce = (el-1)*discretization_nIPs + ip
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
en = material_phaseEntry(co,ce)

6
src/phase_mechanical.f90
View File

@ -84,7 +84,7 @@ submodule(phase) mechanical
ph, &
en
real(pReal), intent(in) :: &
subdt !< timestep
subdt !< timestep
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function plastic_dotState
@ -848,10 +848,10 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ dotState * Delta_t
broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),ph,en)
broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage), ph,en)
if(broken) exit
dotState = plastic_dotState(Delta_t,ph,en)
dotState = plastic_dotState(Delta_t*C(stage), ph,en)
if (any(IEEE_is_NaN(dotState))) exit
enddo

2
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -27,7 +27,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstances,p,i,k
integer :: Ninstances, p, k
class(tNode), pointer :: &
phases, &
phase, &

2
src/phase_mechanical_elastic.f90
View File

@ -48,7 +48,7 @@ module subroutine elastic_init(phases)
prm%C_11 = polynomial(elastic%asDict(),'C_11','T')
prm%C_12 = polynomial(elastic%asDict(),'C_12','T')
prm%C_44 = polynomial(elastic%asDict(),'C_44','T')
if (any(phase_lattice(ph) == ['hP','tI'])) then
prm%C_13 = polynomial(elastic%asDict(),'C_13','T')
prm%C_33 = polynomial(elastic%asDict(),'C_33','T')

44
src/phase_mechanical_plastic.f90
View File

@ -110,29 +110,35 @@ submodule(phase:mechanical) plastic
end subroutine nonlocal_LpAndItsTangent
module subroutine isotropic_dotState(Mp,ph,en)
module function isotropic_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine isotropic_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function isotropic_dotState
module subroutine phenopowerlaw_dotState(Mp,ph,en)
module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine phenopowerlaw_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function phenopowerlaw_dotState
module subroutine plastic_kinehardening_dotState(Mp,ph,en)
module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
end subroutine plastic_kinehardening_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function plastic_kinehardening_dotState
module subroutine dislotwin_dotState(Mp,T,ph,en)
real(pReal), dimension(3,3), intent(in) :: &
@ -144,15 +150,15 @@ submodule(phase:mechanical) plastic
en
end subroutine dislotwin_dotState
module subroutine dislotungsten_dotState(Mp,T,ph,en)
module function dislotungsten_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
ph, &
en
end subroutine dislotungsten_dotState
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function dislotungsten_dotState
module subroutine nonlocal_dotState(Mp,timestep,ph,en)
real(pReal), dimension(3,3), intent(in) :: &
@ -311,27 +317,28 @@ module function plastic_dotState(subdt,ph,en) result(dotState)
plasticType: select case (phase_plasticity(ph))
case (PLASTIC_ISOTROPIC_ID) plasticType
call isotropic_dotState(Mp,ph,en)
dotState = isotropic_dotState(Mp,ph,en)
case (PLASTIC_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_dotState(Mp,ph,en)
dotState = phenopowerlaw_dotState(Mp,ph,en)
case (PLASTIC_KINEHARDENING_ID) plasticType
call plastic_kinehardening_dotState(Mp,ph,en)
dotState = plastic_kinehardening_dotState(Mp,ph,en)
case (PLASTIC_DISLOTWIN_ID) plasticType
call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en)
dotState = plasticState(ph)%dotState(:,en)
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en)
dotState = dislotungsten_dotState(Mp,ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType
call nonlocal_dotState(Mp,subdt,ph,en)
dotState = plasticState(ph)%dotState(:,en)
end select plasticType
end if
dotState = plasticState(ph)%dotState(:,en)
end function plastic_dotState
@ -400,10 +407,9 @@ module function plastic_deltaState(ph, en) result(broken)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
if (.not. broken) then
myOffset = plasticState(ph)%offsetDeltaState
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = &
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en)
plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) &
+ plasticState(ph)%deltaState(1:mySize,en)
end if
end select

57
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -43,6 +43,13 @@ submodule(phase:plastic) dislotungsten
systems_sl
end type tParameters !< container type for internal constitutive parameters
type :: tIndexDotState
integer, dimension(2) :: &
rho_mob, &
rho_dip, &
gamma_sl
end type tIndexDotState
type :: tDislotungstenState
real(pReal), dimension(:,:), pointer :: &
rho_mob, &
@ -58,10 +65,9 @@ submodule(phase:plastic) dislotungsten
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tDisloTungstenState), allocatable, dimension(:) :: &
dotState, &
state
type(tParameters), allocatable, dimension(:) :: param
type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tDisloTungstenState), allocatable, dimension(:) :: state
type(tDisloTungstenDependentState), allocatable, dimension(:) :: dependentState
contains
@ -103,18 +109,17 @@ module function plastic_dislotungsten_init() result(myPlasticity)
print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002'
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(dependentState(phases%length))
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
@ -214,28 +219,29 @@ module function plastic_dislotungsten_init() result(myPlasticity)
sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
idx_dot%rho_mob = [startIndex,endIndex]
stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob = spread(rho_mob_0,2,Nmembers)
dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
idx_dot%rho_dip = [startIndex,endIndex]
stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip = spread(rho_dip_0,2,Nmembers)
dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
idx_dot%gamma_sl = [startIndex,endIndex]
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
@ -300,15 +306,15 @@ end subroutine dislotungsten_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine dislotungsten_dotState(Mp,T,ph,en)
module function dislotungsten_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T !< temperature
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_pos, dot_gamma_neg,&
@ -319,17 +325,22 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot_rho_dip_climb, &
d_hat
real(pReal) :: &
mu
mu, T
associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
dot_rho_mob => dotState(indexDotState(ph)%rho_mob(1):indexDotState(ph)%rho_mob(2)), &
dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), &
dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)))
mu = elastic_mu(ph,en)
T = thermal_T(ph,en)
call kinetics(Mp,T,ph,en,&
dot_gamma_pos,dot_gamma_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg)
dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
dot_gamma_sl = abs(dot_gamma_pos+dot_gamma_neg)
where(dEq0((tau_pos+tau_neg)*0.5_pReal))
dot_rho_dip_formation = 0.0_pReal
@ -338,7 +349,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
d_hat = math_clip(3.0_pReal*mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), &
prm%d_caron, & ! lower limit
dst%Lambda_sl(:,en)) ! upper limit
dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, &
dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot_gamma_sl/prm%b_sl, &
0.0_pReal, &
prm%dipoleformation)
v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(TAU*K_B*T)) &
@ -346,16 +357,16 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency?
end where
dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication
dot_rho_mob = dot_gamma_sl/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication
- dot_rho_dip_formation &
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges
dot%rho_dip(:,en) = dot_rho_dip_formation &
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot_gamma_sl ! Spontaneous annihilation of 2 edges
dot_rho_dip = dot_rho_dip_formation &
- (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot_gamma_sl & ! Spontaneous annihilation of an edge with a dipole
- dot_rho_dip_climb
end associate
end subroutine dislotungsten_dotState
end function dislotungsten_dotState
!--------------------------------------------------------------------------------------------------

44
src/phase_mechanical_plastic_isotropic.f90
View File

@ -37,9 +37,7 @@ submodule(phase:plastic) isotropic
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tIsotropicState), allocatable, dimension(:) :: &
dotState, &
state
type(tIsotropicState), allocatable, dimension(:) :: state
contains
@ -77,16 +75,15 @@ module function plastic_isotropic_init() result(myPlasticity)
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
associate(prm => param(ph), stt => state(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
#if defined (__GFORTRAN__)
prm%output = output_as1dString(pl)
@ -125,12 +122,12 @@ module function plastic_isotropic_init() result(myPlasticity)
sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
stt%xi => plasticState(ph)%state (1,:)
stt%xi = xi_0
dot%xi => plasticState(ph)%dotState(1,:)
stt%xi => plasticState(ph)%state(1,:)
stt%xi = xi_0
plasticState(ph)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if (plasticState(ph)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
@ -178,7 +175,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
if (norm_Mp_dev > 0.0_pReal) then
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en)))**prm%n
Lp = dot_gamma * Mp_dev/norm_Mp_dev
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -242,27 +239,26 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine isotropic_dotState(Mp,ph,en)
module function isotropic_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal) :: &
dot_gamma, & !< strainrate
xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(ph), stt => state(ph), &
dot => dotState(ph))
associate(prm => param(ph), stt => state(ph), dot_xi => dotState(1))
if (prm%dilatation) then
norm_Mp = sqrt(math_tensordot(Mp,Mp))
else
norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp)))
end if
norm_Mp = merge(sqrt(math_tensordot(Mp,Mp)), &
sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))), &
prm%dilatation)
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n
@ -274,16 +270,16 @@ module subroutine isotropic_dotState(Mp,ph,en)
+ asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) &
/ prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
end if
dot%xi(en) = dot_gamma &
* ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star)
dot_xi = dot_gamma &
* ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star)
else
dot%xi(en) = 0.0_pReal
dot_xi = 0.0_pReal
end if
end associate
end subroutine isotropic_dotState
end function isotropic_dotState
!--------------------------------------------------------------------------------------------------

56
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -34,6 +34,13 @@ submodule(phase:plastic) kinehardening
systems_sl
end type tParameters
type :: tIndexDotState
integer, dimension(2) :: &
xi, &
chi, &
gamma
end type tIndexDotState
type :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: &
xi, & !< resistance against plastic slip
@ -47,10 +54,8 @@ submodule(phase:plastic) kinehardening
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param
type(tKinehardeningState), allocatable, dimension(:) :: &
dotState, &
deltaState, &
state
type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tKinehardeningState), allocatable, dimension(:) :: state, deltaState
contains
@ -91,19 +96,20 @@ module function plastic_kinehardening_init() result(myPlasticity)
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(deltaState(phases%length))
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph))
associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
#if defined (__GFORTRAN__)
prm%output = output_as1dString(pl)
@ -173,27 +179,28 @@ module function plastic_kinehardening_init() result(myPlasticity)
sizeState = sizeDotState + sizeDeltaState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi => plasticState(ph)%state (startIndex:endIndex,:)
idx_dot%xi = [startIndex,endIndex]
stt%xi => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi = spread(xi_0, 2, Nmembers)
dot%xi => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%chi => plasticState(ph)%state (startIndex:endIndex,:)
dot%chi => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%chi = [startIndex,endIndex]
stt%chi => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%gamma = [startIndex,endIndex]
stt%gamma => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
@ -270,13 +277,15 @@ end subroutine kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_kinehardening_dotState(Mp,ph,en)
module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal) :: &
sumGamma
@ -284,29 +293,32 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en)
dot_gamma_pos,dot_gamma_neg
associate(prm => param(ph), stt => state(ph),dot => dotState(ph))
associate(prm => param(ph), stt => state(ph), &
dot_xi => dotState(IndexDotState(ph)%xi(1):IndexDotState(ph)%xi(2)),&
dot_chi => dotState(IndexDotState(ph)%chi(1):IndexDotState(ph)%chi(2)),&
dot_gamma => dotState(IndexDotState(ph)%gamma(1):IndexDotState(ph)%gamma(2)))
call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
dot%gamma(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
dot_gamma = abs(dot_gamma_pos+dot_gamma_neg)
sumGamma = sum(stt%gamma(:,en))
dot%xi(:,en) = matmul(prm%h_sl_sl,dot%gamma(:,en)) &
dot_xi = matmul(prm%h_sl_sl,dot_gamma) &
* ( prm%h_inf_f &
+ (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
* exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
)
dot%chi(:,en) = stt%sgn_gamma(:,en)*dot%gamma(:,en) * &
( prm%h_inf_b + &
(prm%h_0_b - prm%h_inf_b &
dot_chi = stt%sgn_gamma(:,en)*dot_gamma &
* ( prm%h_inf_b &
+ (prm%h_0_b - prm%h_inf_b &
+ prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi_0(:,en))*(stt%gamma(:,en)-stt%gamma_0(:,en))&
) *exp(-(stt%gamma(:,en)-stt%gamma_0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi_0(:,en))) &
)
end associate
end subroutine plastic_kinehardening_dotState
end function plastic_kinehardening_dotState
!--------------------------------------------------------------------------------------------------

4
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -231,8 +231,8 @@ module function plastic_nonlocal_init() result(myPlasticity)
st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
plasticState(ph)%nonlocal = pl%get_asBool('flux',defaultVal=.True.)
#if defined (__GFORTRAN__)

70
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -47,6 +47,14 @@ submodule(phase:plastic) phenopowerlaw
systems_tw
end type tParameters
type :: tIndexDotState
integer, dimension(2) :: &
xi_sl, &
xi_tw, &
gamma_sl, &
gamma_tw
end type tIndexDotState
type :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
xi_sl, &
@ -56,11 +64,10 @@ submodule(phase:plastic) phenopowerlaw
end type tPhenopowerlawState
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
! containers for parameters, dot state index, and state
type(tParameters), allocatable, dimension(:) :: param
type(tPhenopowerlawState), allocatable, dimension(:) :: &
dotState, &
state
type(tIndexDotState), allocatable, dimension(:) :: indexDotState
type(tPhenopowerlawState), allocatable, dimension(:) :: state
contains
@ -101,17 +108,18 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
phases => config_material%get('phase')
allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length))
allocate(dotState(phases%length))
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
associate(prm => param(ph), stt => state(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
mech => phase%get('mechanical')
pl => mech%get('plastic')
!--------------------------------------------------------------------------------------------------
! slip related parameters
@ -224,37 +232,37 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi_sl => plasticState(ph)%state (startIndex:endIndex,:)
idx_dot%xi_sl = [startIndex,endIndex]
stt%xi_sl => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_sl = spread(xi_0_sl, 2, Nmembers)
dot%xi_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%xi_tw => plasticState(ph)%state (startIndex:endIndex,:)
idx_dot%xi_tw = [startIndex,endIndex]
stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_tw = spread(xi_0_tw, 2, Nmembers)
dot%xi_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%gamma_sl = [startIndex,endIndex]
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%gamma_tw => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
idx_dot%gamma_tw = [startIndex,endIndex]
stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
end associate
@ -324,13 +332,15 @@ end subroutine phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
module subroutine phenopowerlaw_dotState(Mp,ph,en)
module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
ph, &
en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal) :: &
xi_sl_sat_offset,&
@ -340,28 +350,32 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en)
right_SlipSlip
associate(prm => param(ph), stt => state(ph), dot => dotState(ph))
associate(prm => param(ph), stt => state(ph), &
dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), &
dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), &
dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), &
dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2)))
call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en))
dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
call kinetics_tw(Mp,ph,en,dot_gamma_tw)
sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, &
1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
dot%xi_sl(:,en) = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
* matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot%gamma_tw(:,en))
dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
* matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot_gamma_tw)
dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
* matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) &
+ prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en))
dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
* matmul(prm%h_tw_sl,dot_gamma_sl) &
+ prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw)
end associate
end subroutine phenopowerlaw_dotState
end function phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------

21
src/phase_thermal.f90
View File

@ -96,17 +96,24 @@ module subroutine thermal_init(phases)
do ph = 1, phases%length
Nmembers = count(material_phaseID == ph)
allocate(current(ph)%T(Nmembers),source=300.0_pReal)
allocate(current(ph)%T(Nmembers),source=T_ROOM)
allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
phase => phases%get(ph)
thermal => phase%get('thermal',defaultVal=emptyDict)
param(ph)%C_p = thermal%get_asFloat('C_p',defaultVal=0.0_pReal) ! ToDo: make mandatory?
param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal) ! ToDo: make mandatory?
param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal) ! ToDo: depends on symmtery
param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length
! ToDo: temperature dependency of K and C_p
if (thermal%length > 0) then
param(ph)%C_p = thermal%get_asFloat('C_p')
param(ph)%K(1,1) = thermal%get_asFloat('K_11')
if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asFloat('K_33')
param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length
else
thermal_Nsources(ph) = 0
end if
allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
enddo

2
src/prec.f90
View File

@ -45,6 +45,8 @@ module prec
state, & !< state
dotState, & !< rate of state change
deltaState !< increment of state change
real(pReal), pointer, dimension(:,:) :: &
deltaState2
end type
type, extends(tState) :: tPlasticState

23
src/system_routines.f90
View File

@ -24,7 +24,7 @@ module system_routines
function setCWD_C(cwd) bind(C)
use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
integer(C_INT) :: setCWD_C
character(kind=C_CHAR), dimension(*), intent(in) :: cwd
end function setCWD_C
@ -150,14 +150,14 @@ function getUserName()
getUserName = c_f_string(getUserName_Cstring)
else
getUserName = 'n/a (Error!)'
endif
end if
end function getUserName
!--------------------------------------------------------------------------------------------------
!> @brief convert C string to Fortran string
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string
!> @brief Convert C string to Fortran string.
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string.
!--------------------------------------------------------------------------------------------------
pure function c_f_string(c_string) result(f_string)
@ -174,28 +174,23 @@ pure function c_f_string(c_string) result(f_string)
else
f_string = f_string(:i-1)
exit
endif
enddo arrayToString
end if
end do arrayToString
end function c_f_string
!--------------------------------------------------------------------------------------------------
!> @brief convert Fortran string to C string
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string
!> @brief Convert Fortran string to C string.
!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string.
!--------------------------------------------------------------------------------------------------
pure function f_c_string(f_string) result(c_string)
character(len=*), intent(in) :: f_string
character(kind=C_CHAR), dimension(len_trim(f_string)+1) :: c_string
integer :: i
do i=1,len_trim(f_string)
c_string(i)=f_string(i:i)
enddo
c_string(len_trim(f_string)+1) = C_NULL_CHAR
c_string = transfer(trim(f_string)//C_NULL_CHAR,c_string,size=size(c_string))
end function f_c_string