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 untangle-homog-phase

This commit is contained in:
Martin Diehl 2023-07-18 21:15:46 +02:00
parent dc9d4bb5a9 5bd9531577
commit 5e9162b073
33 changed files with 580 additions and 490 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
.gitlab-ci.yml
View File

@ -209,7 +209,7 @@ grid_performance:
- make -j2 all install - make -j2 all install
- export PATH=${PWD}/bin:${PATH} - export PATH=${PWD}/bin:${PATH}
- cd $(mktemp -d) - cd $(mktemp -d)
- git clone -q git@git.damask.mpie.de:damask/performance.git . - git clone -q git@git.damask.mpie.de:damask/statistics.git .
- ./measure_performance.py --input_dir ${CI_PROJECT_DIR}/examples/grid --tag ${CI_COMMIT_SHA} - ./measure_performance.py --input_dir ${CI_PROJECT_DIR}/examples/grid --tag ${CI_COMMIT_SHA}
- > - >
if [ ${CI_COMMIT_BRANCH} == development ]; then if [ ${CI_COMMIT_BRANCH} == development ]; then
@ -222,7 +222,7 @@ update_plots:
stage: statistics stage: statistics
script: script:
- cd $(mktemp -d) - cd $(mktemp -d)
- git clone -q git@git.damask.mpie.de:damask/performance.git . - git clone -q git@git.damask.mpie.de:damask/statistics.git .
- ./plot_commithistory.py --color green -n 5 -N 100 - ./plot_commithistory.py --color green -n 5 -N 100
- ./plot_commithistory.py --color green -n 5 -N 1000 - ./plot_commithistory.py --color green -n 5 -N 1000
- ./plot_commithistory.py --color green -n 5 -N 10000 - ./plot_commithistory.py --color green -n 5 -N 10000

2
PRIVATE

@ -1 +1 @@
Subproject commit 39853f7d8fb06bb0c83aed5083212bcfd5b52a57 Subproject commit 162106e6379d484ee101981c66e3f159d2f8821a

2
VERSION
View File

@ -1 +1 @@
3.0.0-alpha7-630-ga63fe8d02 3.0.0-alpha7-665-g341119d70

50
examples/config/numerics.yaml
View File

@ -25,30 +25,40 @@ homogenization:
stepIncrease: 1.5 # increase of next substep size when previous substep converged in homogenization (value higher than 1) stepIncrease: 1.5 # increase of next substep size when previous substep converged in homogenization (value higher than 1)
nMPstate: 10 # materialpoint state loop limit nMPstate: 10 # materialpoint state loop limit
grid: solver:
eps_div_atol: 1.0e-3 # absolute tolerance for fulfillment of stress equilibrium grid:
eps_div_rtol: 5.0e-4 # relative tolerance for fulfillment of stress equilibrium N_staggered_iter_max: 10 # max number of field level staggered iterations
eps_curl_atol: 1.0e-12 # absolute tolerance for fulfillment of strain compatibility N_cutback_max: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
eps_curl_rtol: 5.0e-4 # relative tolerance for fulfillment of strain compatibility
eps_stress_atol: 1.0e+3 # absolute tolerance for fulfillment of stress BC damage:
eps_stress_rtol: 0.01 # relative tolerance for fulfillment of stress BC N_iter_max: 100 # maximum iteration number
eps_damage_atol: 1.0e-2 # absolute tolerance for damage evolution eps_abs_phi: 1.0e-2 # absolute tolerance for damage evolution
eps_damage_rtol: 1.0e-6 # relative tolerance for damage evolution eps_rel_phi: 1.0e-6 # relative tolerance for damage evolution
eps_thermal_atol: 1.0e-2 # absolute tolerance for thermal equilibrium thermal:
eps_thermal_rtol: 1.0e-6 # relative tolerance for thermal equilibrium N_iter_max: 100 # maximum iteration number
itmax: 250 # Maximum iteration number eps_abs_T: 1.0e-2 # absolute tolerance for thermal equilibrium
itmin: 2 # Minimum iteration number eps_rel_T: 1.0e-6 # relative tolerance for thermal equilibrium
fftw_timelimit: -1.0 # timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit
fftw_plan_mode: FFTW_PATIENT # reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag mechanical:
maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc) eps_abs_div(P): 1.0e-4 # absolute tolerance for fulfillment of stress equilibrium
maxStaggeredIter: 10 # max number of field level staggered iterations eps_rel_div(P): 5.0e-4 # relative tolerance for fulfillment of stress equilibrium
memory_efficient: 1 # Precalculate Gamma-operator (81 double per point) eps_abs_P: 1.0e3 # absolute tolerance for fulfillment of stress BC
eps_rel_P: 1.0e-3 # relative tolerance for fulfillment of stress BC
N_iter_min: 1 # minimum iteration number
N_iter_max: 100 # maximum iteration number
update_gamma: false # Update Gamma-operator with current dPdF (not possible if memory_efficient=1) update_gamma: false # Update Gamma-operator with current dPdF (not possible if memory_efficient=1)
divergence_correction: 2 # Use size-independent divergence criterion
FFT:
memory_efficient: True # Precalculate Gamma-operator (81 double per point)
divergence_correction: size+grid # Use size-independent divergence criterion
derivative: continuous # Approximation used for derivatives in Fourier space derivative: continuous # Approximation used for derivatives in Fourier space
petsc_options: -snes_type ngmres -snes_ngmres_anderson # PetSc solver options FFTW_plan_mode: FFTW_MEASURE # planing-rigor flag, see manual on www.fftw.org
FFTW_timelimit: -1.0 # timelimit of plan creation for FFTW, see manual on www.fftw.org. -1.0: disable timelimit
PETSc_options: -snes_type ngmres -snes_ngmres_anderson # PETSc solver options
alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
eps_abs_curl(F): 1.0e-10 # absolute tolerance for fulfillment of strain compatibility
eps_rel_curl(F): 5.0e-4 # relative tolerance for fulfillment of strain compatibility
mesh: mesh:
maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc) maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)

8
examples/grid/numerics.yaml
View File

@ -1,3 +1,5 @@
grid: solver:
itmin: 4 grid:
itmax: 40 mechanical:
N_iter_min: 4
N_iter_max: 40

57
src/CLI.f90
View File

@ -74,7 +74,7 @@ subroutine CLI_init()
print'(a)', achar(27)//'[31m' print'(a)', achar(27)//'[31m'
print'(1x,a,/)', 'debug version - debug version - debug version - debug version - debug version' print'(1x,a,/)', 'debug version - debug version - debug version - debug version - debug version'
#else #else
print'(a)', achar(27)//'[94m' print'(a)', achar(27)//'[1;94m'
#endif #endif
print'(1x,a)', ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/ _/_/_/' print'(1x,a)', ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/ _/_/_/'
print'(1x,a)', ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/ _/' print'(1x,a)', ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/ _/'
@ -158,48 +158,34 @@ subroutine CLI_init()
print'(1x,a,/)',' Prints this message and exits' print'(1x,a,/)',' Prints this message and exits'
call quit(0) ! normal Termination call quit(0) ! normal Termination
case ('-g', '--geom', '--geometry') case ('-g', '--geom', '--geometry')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --geom' if (.not. hasArg) call IO_error(610,ext_msg='--geom')
geomArg = getArg(i+1) geomArg = getArg(i+1)
case ('-l', '--load', '--loadcase') case ('-l', '--load', '--loadcase')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --load' if (.not. hasArg) call IO_error(610,ext_msg='--load')
loadArg = getArg(i+1) loadArg = getArg(i+1)
case ('-m', '--material', '--materialconfig') case ('-m', '--material', '--materialconfig')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --material' if (.not. hasArg) call IO_error(610,ext_msg='--material')
materialArg = getArg(i+1) materialArg = getArg(i+1)
case ('-n', '--numerics', '--numericsconfig') case ('-n', '--numerics', '--numericsconfig')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --numerics' if (.not. hasArg) call IO_error(610,ext_msg='--numerics')
numericsArg = getArg(i+1) numericsArg = getArg(i+1)
case ('-j', '--job', '--jobname') case ('-j', '--job', '--jobname')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --jobname' if (.not. hasArg) call IO_error(610,ext_msg='--jobname')
solverJobname = getArg(i+1) solverJobname = getArg(i+1)
case ('-w', '--wd', '--workingdir', '--workingdirectory') case ('-w', '--wd', '--workingdir', '--workingdirectory')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --workingdirectory' if (.not. hasArg) call IO_error(610,ext_msg='--workingdirectory')
workingDirArg = getArg(i+1) workingDirArg = getArg(i+1)
case ('-r', '--rs', '--restart') case ('-r', '--rs', '--restart')
if (.not. hasArg) print'(/,1x,a)', 'ERROR: missing argument for --restart' if (.not. hasArg) call IO_error(610,ext_msg='--jobname')
arg = getArg(i+1) arg = getArg(i+1)
read(arg,*,iostat=stat) CLI_restartInc read(arg,*,iostat=stat) CLI_restartInc
if (CLI_restartInc < 0 .or. stat /= 0) then if (CLI_restartInc < 0 .or. stat /= 0) call IO_error(611,ext_msg=arg)
print'(/,1x,a)', 'ERROR: could not parse restart increment: '//trim(arg)
call quit(1)
end if
end select end select
end do end do
if (.not. allocated(loadArg)) then if (.not. allocated(geomArg)) call IO_error(612,ext_msg='--geom')
print'(/,1x,a)', 'Error: no load case specified (-h for help)' if (.not. allocated(loadArg)) call IO_error(612,ext_msg='--load')
call quit(1) if (.not. allocated(materialArg)) call IO_error(612,ext_msg='--material')
end if
if (.not. allocated(geomArg)) then
print'(/,1x,a)', 'Error: no geometry specified (-h for help)'
call quit(1)
end if
if (.not. allocated(materialArg)) then
print'(/,1x,a)', 'Error: no material configuration specified (-h for help)'
call quit(1)
end if
call setWorkingDirectory(trim(workingDirArg)) call setWorkingDirectory(trim(workingDirArg))
CLI_geomFile = getPathRelCWD(geomArg,'geometry') CLI_geomFile = getPathRelCWD(geomArg,'geometry')
@ -211,8 +197,7 @@ subroutine CLI_init()
if (.not. allocated(solverJobname)) then if (.not. allocated(solverJobname)) then
solverJobname = jobname(CLI_geomFile,CLI_loadFile,CLI_materialFile,CLI_numericsFile) solverJobname = jobname(CLI_geomFile,CLI_loadFile,CLI_materialFile,CLI_numericsFile)
elseif (scan(solverJobname,'/') > 0) then elseif (scan(solverJobname,'/') > 0) then
print'(/,1x,a)', 'ERROR: JOBNAME must not contain any slashes' call IO_error(630)
call quit(1)
endif endif
commandLine = getArg(-1) commandLine = getArg(-1)
@ -272,9 +257,6 @@ subroutine setWorkingDirectory(workingDirectoryArg)
character(len=*), intent(in) :: workingDirectoryArg !< working directory argument character(len=*), intent(in) :: workingDirectoryArg !< working directory argument
character(len=:), allocatable :: workingDirectory character(len=:), allocatable :: workingDirectory
logical :: error
external :: quit
absolutePath: if (workingDirectoryArg(1:1) == '/') then absolutePath: if (workingDirectoryArg(1:1) == '/') then
workingDirectory = workingDirectoryArg workingDirectory = workingDirectoryArg
@ -284,11 +266,7 @@ subroutine setWorkingDirectory(workingDirectoryArg)
end if absolutePath end if absolutePath
workingDirectory = trim(normpath(workingDirectory)) workingDirectory = trim(normpath(workingDirectory))
error = setCWD(trim(workingDirectory)) if (setCWD(trim(workingDirectory))) call IO_error(640,ext_msg=workingDirectory)
if (error) then
print'(1x,a)', 'ERROR: invalid working directory: '//trim(workingDirectory)
call quit(1)
end if
end subroutine setWorkingDirectory end subroutine setWorkingDirectory
@ -344,7 +322,6 @@ function getPathRelCWD(path,fileType)
character(len=*), intent(in) :: fileType character(len=*), intent(in) :: fileType
logical :: file_exists logical :: file_exists
external :: quit
getPathRelCWD = trim(path) getPathRelCWD = trim(path)
@ -352,10 +329,7 @@ function getPathRelCWD(path,fileType)
getPathRelCWD = trim(relpath(getPathRelCWD,getCWD())) getPathRelCWD = trim(relpath(getPathRelCWD,getCWD()))
inquire(file=getPathRelCWD, exist=file_exists) inquire(file=getPathRelCWD, exist=file_exists)
if (.not. file_exists) then if (.not. file_exists) call IO_error(100,ext_msg=fileType//' "'//trim(getPathRelCWD)//'"')
print'(/,1x,a)', 'ERROR: '//fileType//' file does not exist: '//trim(getPathRelCWD)
call quit(1)
end if
end function getPathRelCWD end function getPathRelCWD
@ -438,4 +412,5 @@ function relpath(path,start)
end function relpath end function relpath
end module CLI end module CLI

2
src/CMakeLists.txt
View File

@ -1,7 +1,7 @@
# special flags for some files # special flags for some files
if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
# long lines for interaction matrix # long lines for interaction matrix
set_source_files_properties("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240") set_source_files_properties("crystal.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240")
set_source_files_properties("parallelization.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-none") set_source_files_properties("parallelization.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-none")
endif() endif()

56
src/IO.f90
View File

@ -11,6 +11,7 @@ module IO
IO_STDERR => ERROR_UNIT IO_STDERR => ERROR_UNIT
use prec use prec
use constants
use misc use misc
implicit none(type,external) implicit none(type,external)
@ -20,11 +21,8 @@ module IO
IO_WHITESPACE = achar(44)//achar(32)//achar(9)//achar(10)//achar(13), & !< whitespace characters IO_WHITESPACE = achar(44)//achar(32)//achar(9)//achar(10)//achar(13), & !< whitespace characters
IO_QUOTES = "'"//'"' IO_QUOTES = "'"//'"'
character, parameter, public :: & character, parameter, public :: &
IO_EOL = new_line('DAMASK'), & !< end of line character IO_EOL = LF, & !< end of line character
IO_COMMENT = '#' IO_COMMENT = '#'
character, parameter :: &
CR = achar(13), &
LF = IO_EOL
public :: & public :: &
IO_init, & IO_init, &
@ -294,9 +292,6 @@ pure function IO_lc(str)
character(len=*), intent(in) :: str !< string to convert character(len=*), intent(in) :: str !< string to convert
character(len=len(str)) :: IO_lc character(len=len(str)) :: IO_lc
character(len=*), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz'
character(len=len(LOWER)), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
integer :: i,n integer :: i,n
@ -476,7 +471,7 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2)
case (131) case (131)
msg = 'hex lattice structure with invalid c/a ratio' msg = 'hex lattice structure with invalid c/a ratio'
case (132) case (132)
msg = 'trans_lattice_structure not possible' msg = 'invalid parameters for transformation'
case (134) case (134)
msg = 'negative lattice parameter' msg = 'negative lattice parameter'
case (135) case (135)
@ -555,6 +550,18 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2)
! user errors ! user errors
case (603) case (603)
msg = 'invalid data for table' msg = 'invalid data for table'
case (610)
msg = 'missing argument for option'
case (611)
msg = 'could not parse restart increment'
case (612)
msg = 'missing option'
case (630)
msg = 'JOBNAME must not contain any slashes'
case (640)
msg = 'invalid working directory'
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
! errors related to YAML data ! errors related to YAML data
@ -704,38 +711,43 @@ subroutine panel(paneltype,ID,msg,ext_msg,label1,ID1,label2,ID2)
character(len=pSTRLEN) :: formatString character(len=pSTRLEN) :: formatString
integer, parameter :: panelwidth = 69 integer, parameter :: panelwidth = 69
character(len=:), allocatable :: msg_,ID_,msg1,msg2
character(len=*), parameter :: DIVIDER = repeat('─',panelwidth) character(len=*), parameter :: DIVIDER = repeat('─',panelwidth)
if (.not. present(label1) .and. present(ID1)) error stop 'missing label for value 1' if (.not. present(label1) .and. present(ID1)) error stop 'missing label for value 1'
if (.not. present(label2) .and. present(ID2)) error stop 'missing label for value 2' if (.not. present(label2) .and. present(ID2)) error stop 'missing label for value 2'
if ( present(label1) .and. .not. present(ID1)) error stop 'missing value for label 1'
if ( present(label2) .and. .not. present(ID2)) error stop 'missing value for label 2'
ID_ = IO_intAsStr(ID)
if (present(label1)) msg1 = label1
if (present(label2)) msg2 = label2
if (present(ID1)) msg1 = msg1//' '//IO_intAsStr(ID1)
if (present(ID2)) msg2 = msg2//' '//IO_intAsStr(ID2)
if (paneltype == 'error') msg_ = achar(27)//'[31m'//trim(msg)//achar(27)//'[0m'
if (paneltype == 'warning') msg_ = achar(27)//'[33m'//trim(msg)//achar(27)//'[0m'
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(IO_STDERR,'(/,a)') ' ┌'//DIVIDER//'┐' write(IO_STDERR,'(/,a)') ' ┌'//DIVIDER//'┐'
write(formatString,'(a,i2,a)') '(a,24x,a,',max(1,panelwidth-24-len_trim(paneltype)),'x,a)' write(formatString,'(a,i2,a)') '(a,24x,a,1x,i0,',max(1,panelwidth-24-len_trim(paneltype)-1-len_trim(ID_)),'x,a)'
write(IO_STDERR,formatString) ' │',trim(paneltype), '│' write(IO_STDERR,formatString) ' │',trim(paneltype),ID, '│'
write(formatString,'(a,i2,a)') '(a,24x,i3,',max(1,panelwidth-24-3),'x,a)'
write(IO_STDERR,formatString) ' │',ID, '│'
write(IO_STDERR,'(a)') ' ├'//DIVIDER//'┤' write(IO_STDERR,'(a)') ' ├'//DIVIDER//'┤'
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(msg)),',',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(msg_)),',',&
max(1,panelwidth+3-len_trim(msg)-4),'x,a)' max(1,panelwidth+3-len_trim(msg)-4),'x,a)'
write(IO_STDERR,formatString) '│ ',trim(msg), '│' write(IO_STDERR,formatString) '│ ',trim(msg_), '│'
if (present(ext_msg)) then if (present(ext_msg)) then
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
max(1,panelwidth+3-len_trim(ext_msg)-4),'x,a)' max(1,panelwidth+3-len_trim(ext_msg)-4),'x,a)'
write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│' write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│'
end if end if
if (present(label1)) then if (present(label1)) then
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label1)),',i9,',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(msg1)),',',&
max(1,panelwidth+3-len_trim(label1)-9-7),'x,a)' max(1,panelwidth+3-len_trim(msg1)-7),'x,a)'
write(IO_STDERR,formatString) '│ at ',trim(label1),ID1, '│' write(IO_STDERR,formatString) '│ at ',trim(msg1), '│'
end if end if
if (present(label2)) then if (present(label2)) then
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label2)),',i9,',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(msg2)),',',&
max(1,panelwidth+3-len_trim(label2)-9-7),'x,a)' max(1,panelwidth+3-len_trim(msg2)-7),'x,a)'
write(IO_STDERR,formatString) '│ at ',trim(label2),ID2, '│' write(IO_STDERR,formatString) '│ at ',trim(msg2), '│'
end if end if
write(formatString,'(a,i2.2,a)') '(a,',max(1,panelwidth),'x,a)' write(formatString,'(a,i2.2,a)') '(a,',max(1,panelwidth),'x,a)'
write(IO_STDERR,formatString) ' │', '│' write(IO_STDERR,formatString) ' │', '│'

4
src/Marc/DAMASK_Marc.f90
View File

@ -142,8 +142,8 @@ end function solverIsSymmetric
end module DAMASK_interface end module DAMASK_interface
#include "../parallelization.f90" #include "../parallelization.f90"
#include "../misc.f90"
#include "../constants.f90" #include "../constants.f90"
#include "../misc.f90"
#include "../IO.f90" #include "../IO.f90"
#include "../YAML_types.f90" #include "../YAML_types.f90"
#include "../YAML_parse.f90" #include "../YAML_parse.f90"
@ -155,7 +155,7 @@ end module DAMASK_interface
#include "../rotations.f90" #include "../rotations.f90"
#include "../polynomials.f90" #include "../polynomials.f90"
#include "../tables.f90" #include "../tables.f90"
#include "../lattice.f90" #include "../crystal.f90"
#include "element.f90" #include "element.f90"
#include "../geometry_plastic_nonlocal.f90" #include "../geometry_plastic_nonlocal.f90"
#include "../discretization.f90" #include "../discretization.f90"

4
src/Marc/materialpoint_Marc.f90
View File

@ -16,7 +16,7 @@ module materialpoint_Marc
use rotations use rotations
use polynomials use polynomials
use tables use tables
use lattice use crystal
use material use material
use phase use phase
use homogenization use homogenization
@ -75,7 +75,7 @@ subroutine materialpoint_initAll()
call rotations_init() call rotations_init()
call polynomials_init() call polynomials_init()
call tables_init() call tables_init()
call lattice_init() call crystal_init()
call discretization_Marc_init() call discretization_Marc_init()
call material_init(.false.) call material_init(.false.)
call phase_init() call phase_init()

7
src/constants.f90
View File

@ -13,4 +13,11 @@ module constants
K_B = 1.380649e-23_pREAL, & !< Boltzmann constant in J/Kelvin (https://doi.org/10.1351/goldbook) K_B = 1.380649e-23_pREAL, & !< Boltzmann constant in J/Kelvin (https://doi.org/10.1351/goldbook)
N_A = 6.02214076e23_pREAL !< Avogadro constant in 1/mol (https://doi.org/10.1351/goldbook) N_A = 6.02214076e23_pREAL !< Avogadro constant in 1/mol (https://doi.org/10.1351/goldbook)
character, parameter :: &
CR = achar(13), &
LF = new_line('DAMASK')
character(len=*), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz'
character(len=len(LOWER)), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
end module constants end module constants

316
src/lattice.f90 → src/crystal.f90
View File

@ -3,10 +3,10 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief contains lattice definitions including Schmid matrices for slip, twin, trans, !> @brief Contains crystal definitions including Schmid matrices for slip, twin, trans,
! and cleavage as well as interaction among the various systems ! and cleavage as well as interaction among the various systems.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module lattice module crystal
use prec use prec
use misc use misc
use IO use IO
@ -80,7 +80,7 @@ module lattice
],pREAL),shape(CF_SYSTEMTWIN)) !< cF twin systems ],pREAL),shape(CF_SYSTEMTWIN)) !< cF twin systems
integer, dimension(2,CF_NTWIN), parameter, public :: & integer, dimension(2,CF_NTWIN), parameter, public :: &
lattice_CF_TWINNUCLEATIONSLIPPAIR = reshape( [& crystal_CF_TWINNUCLEATIONSLIPPAIR = reshape( [&
2,3, & 2,3, &
1,3, & 1,3, &
1,2, & 1,2, &
@ -93,7 +93,7 @@ module lattice
11,12, & 11,12, &
10,12, & 10,12, &
10,11 & 10,11 &
],shape(lattice_CF_TWINNUCLEATIONSLIPPAIR)) ],shape(crystal_CF_TWINNUCLEATIONSLIPPAIR))
real(pREAL), dimension(3+3,CF_NCLEAVAGE), parameter :: & real(pREAL), dimension(3+3,CF_NCLEAVAGE), parameter :: &
CF_SYSTEMCLEAVAGE = reshape(real([& CF_SYSTEMCLEAVAGE = reshape(real([&
@ -367,60 +367,60 @@ module lattice
],pREAL),shape(TI_SYSTEMSLIP)) !< tI slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x) ],pREAL),shape(TI_SYSTEMSLIP)) !< tI slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x)
interface lattice_forestProjection_edge interface crystal_forestProjection_edge
module procedure slipProjection_transverse module procedure slipProjection_transverse
end interface lattice_forestProjection_edge end interface crystal_forestProjection_edge
interface lattice_forestProjection_screw interface crystal_forestProjection_screw
module procedure slipProjection_direction module procedure slipProjection_direction
end interface lattice_forestProjection_screw end interface crystal_forestProjection_screw
public :: & public :: &
lattice_init, & crystal_init, &
lattice_isotropic_nu, & crystal_isotropic_nu, &
lattice_isotropic_mu, & crystal_isotropic_mu, &
lattice_symmetrize_33, & crystal_symmetrize_33, &
lattice_symmetrize_C66, & crystal_symmetrize_C66, &
lattice_SchmidMatrix_slip, & crystal_SchmidMatrix_slip, &
lattice_SchmidMatrix_twin, & crystal_SchmidMatrix_twin, &
lattice_SchmidMatrix_trans, & crystal_SchmidMatrix_trans, &
lattice_SchmidMatrix_cleavage, & crystal_SchmidMatrix_cleavage, &
lattice_nonSchmidMatrix, & crystal_nonSchmidMatrix, &
lattice_interaction_SlipBySlip, & crystal_interaction_SlipBySlip, &
lattice_interaction_TwinByTwin, & crystal_interaction_TwinByTwin, &
lattice_interaction_TransByTrans, & crystal_interaction_TransByTrans, &
lattice_interaction_SlipByTwin, & crystal_interaction_SlipByTwin, &
lattice_interaction_SlipByTrans, & crystal_interaction_SlipByTrans, &
lattice_interaction_TwinBySlip, & crystal_interaction_TwinBySlip, &
lattice_characteristicShear_Twin, & crystal_characteristicShear_Twin, &
lattice_C66_twin, & crystal_C66_twin, &
lattice_C66_trans, & crystal_C66_trans, &
lattice_forestProjection_edge, & crystal_forestProjection_edge, &
lattice_forestProjection_screw, & crystal_forestProjection_screw, &
lattice_slip_normal, & crystal_slip_normal, &
lattice_slip_direction, & crystal_slip_direction, &
lattice_slip_transverse, & crystal_slip_transverse, &
lattice_labels_slip, & crystal_labels_slip, &
lattice_labels_twin crystal_labels_twin
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Run self test. !> @brief Run self test.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine lattice_init() subroutine crystal_init()
print'(/,1x,a)', '<<<+- lattice init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- crystal init -+>>>'; flush(IO_STDOUT)
call selfTest() call selfTest()
end subroutine lattice_init end subroutine crystal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Characteristic shear for twinning !> @brief Characteristic shear for twinning
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear) function crystal_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -470,7 +470,7 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
characteristicShear(a) = 0.5_pREAL*sqrt(2.0_pREAL) characteristicShear(a) = 0.5_pREAL*sqrt(2.0_pREAL)
case('hP') case('hP')
if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) & if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) &
call IO_error(131,ext_msg='lattice_characteristicShear_Twin') call IO_error(131,ext_msg='crystal_characteristicShear_Twin')
p = sum(HP_NTWINSYSTEM(1:f-1))+s p = sum(HP_NTWINSYSTEM(1:f-1))+s
select case(HP_SHEARTWIN(p)) ! from Christian & Mahajan 1995 p.29 select case(HP_SHEARTWIN(p)) ! from Christian & Mahajan 1995 p.29
case (1) ! <-10.1>{10.2} case (1) ! <-10.1>{10.2}
@ -483,24 +483,24 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
characteristicShear(a) = 2.0_pREAL*(cOverA**2-2.0_pREAL)/3.0_pREAL/cOverA characteristicShear(a) = 2.0_pREAL*(cOverA**2-2.0_pREAL)/3.0_pREAL/cOverA
end select end select
case default case default
call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_characteristicShear_Twin: '//trim(lattice))
end select end select
end do mySystems end do mySystems
end do myFamilies end do myFamilies
end function lattice_characteristicShear_Twin end function crystal_characteristicShear_Twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Rotated elasticity matrices for twinning in 6x6-matrix notation !> @brief Rotated elasticity matrices for twinning in 6x6-matrix notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_twin(Ntwin,C66,lattice,CoverA) function crystal_C66_twin(Ntwin,C66,lattice,CoverA)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
real(pREAL), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix real(pREAL), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix
real(pREAL), intent(in) :: cOverA !< c/a ratio real(pREAL), intent(in) :: cOverA !< c/a ratio
real(pREAL), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin real(pREAL), dimension(6,6,sum(Ntwin)) :: crystal_C66_twin
real(pREAL), dimension(3,3,sum(Ntwin)):: coordinateSystem real(pREAL), dimension(3,3,sum(Ntwin)):: coordinateSystem
type(tRotation) :: R type(tRotation) :: R
@ -518,28 +518,28 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
coordinateSystem = buildCoordinateSystem(Ntwin,HP_NSLIPSYSTEM,HP_SYSTEMTWIN,& coordinateSystem = buildCoordinateSystem(Ntwin,HP_NSLIPSYSTEM,HP_SYSTEMTWIN,&
lattice,cOverA) lattice,cOverA)
case default case default
call IO_error(137,ext_msg='lattice_C66_twin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_C66_twin: '//trim(lattice))
end select end select
do i = 1, sum(Ntwin) do i = 1, sum(Ntwin)
call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg? call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg?
lattice_C66_twin(1:6,1:6,i) = R%rotStiffness(C66) crystal_C66_twin(1:6,1:6,i) = R%rotStiffness(C66)
end do end do
end function lattice_C66_twin end function crystal_C66_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation !> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & function crystal_C66_trans(Ntrans,C_parent66,crystal_target, &
cOverA_trans,a_cF,a_cI) cOverA_trans,a_cF,a_cI)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=*), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: crystal_target !< Bravais lattice (Pearson symbol)
real(pREAL), dimension(6,6), intent(in) :: C_parent66 real(pREAL), dimension(6,6), intent(in) :: C_parent66
real(pREAL), optional, intent(in) :: cOverA_trans, a_cF, a_cI real(pREAL), optional, intent(in) :: cOverA_trans, a_cF, a_cI
real(pREAL), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans real(pREAL), dimension(6,6,sum(Ntrans)) :: crystal_C66_trans
real(pREAL), dimension(6,6) :: C_bar66, C_target_unrotated66 real(pREAL), dimension(6,6) :: C_bar66, C_target_unrotated66
real(pREAL), dimension(3,3,sum(Ntrans)) :: Q,S real(pREAL), dimension(3,3,sum(Ntrans)) :: Q,S
@ -548,11 +548,11 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! elasticity matrix of the target phase in cube orientation ! elasticity matrix of the target phase in cube orientation
if (lattice_target == 'hP' .and. present(cOverA_trans)) then if (crystal_target == 'hP' .and. present(cOverA_trans)) then
! https://doi.org/10.1063/1.1663858 eq. (16), eq. (18), eq. (19) ! https://doi.org/10.1063/1.1663858 eq. (16), eq. (18), eq. (19)
! https://doi.org/10.1016/j.actamat.2016.07.032 eq. (47), eq. (48) ! https://doi.org/10.1016/j.actamat.2016.07.032 eq. (47), eq. (48)
if (cOverA_trans < 1.0_pREAL .or. cOverA_trans > 2.0_pREAL) & if (cOverA_trans < 1.0_pREAL .or. cOverA_trans > 2.0_pREAL) &
call IO_error(131,ext_msg='lattice_C66_trans: '//trim(lattice_target)) call IO_error(131,ext_msg='crystal_C66_trans: '//trim(crystal_target))
C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pREAL*C_parent66(4,4))/2.0_pREAL C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pREAL*C_parent66(4,4))/2.0_pREAL
C_bar66(1,2) = (C_parent66(1,1) + 5.0_pREAL*C_parent66(1,2) - 2.0_pREAL*C_parent66(4,4))/6.0_pREAL C_bar66(1,2) = (C_parent66(1,1) + 5.0_pREAL*C_parent66(1,2) - 2.0_pREAL*C_parent66(4,4))/6.0_pREAL
C_bar66(3,3) = (C_parent66(1,1) + 2.0_pREAL*C_parent66(1,2) + 4.0_pREAL*C_parent66(4,4))/3.0_pREAL C_bar66(3,3) = (C_parent66(1,1) + 2.0_pREAL*C_parent66(1,2) + 4.0_pREAL*C_parent66(4,4))/3.0_pREAL
@ -566,13 +566,13 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
C_target_unrotated66(1,3) = C_bar66(1,3) C_target_unrotated66(1,3) = C_bar66(1,3)
C_target_unrotated66(3,3) = C_bar66(3,3) C_target_unrotated66(3,3) = C_bar66(3,3)
C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2/(0.5_pREAL*(C_bar66(1,1) - C_bar66(1,2))) C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2/(0.5_pREAL*(C_bar66(1,1) - C_bar66(1,2)))
C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP') C_target_unrotated66 = crystal_symmetrize_C66(C_target_unrotated66,'hP')
elseif (lattice_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then elseif (crystal_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then
if (a_cI <= 0.0_pREAL .or. a_cF <= 0.0_pREAL) & if (a_cI <= 0.0_pREAL .or. a_cF <= 0.0_pREAL) &
call IO_error(134,ext_msg='lattice_C66_trans: '//trim(lattice_target)) call IO_error(134,ext_msg='crystal_C66_trans: '//trim(crystal_target))
C_target_unrotated66 = C_parent66 C_target_unrotated66 = C_parent66
else else
call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target)) call IO_error(137,ext_msg='crystal_C66_trans : '//trim(crystal_target))
end if end if
do i = 1,6 do i = 1,6
@ -584,10 +584,10 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
do i = 1,sum(Ntrans) do i = 1,sum(Ntrans)
call R%fromMatrix(Q(1:3,1:3,i)) call R%fromMatrix(Q(1:3,1:3,i))
lattice_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66) crystal_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66)
end do end do
end function lattice_C66_trans end function crystal_C66_trans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -595,7 +595,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
! https://doi.org/10.1016/j.actamat.2012.03.053, eq. (17) ! https://doi.org/10.1016/j.actamat.2012.03.053, eq. (17)
! https://doi.org/10.1016/j.actamat.2008.07.037, table 1 ! https://doi.org/10.1016/j.actamat.2008.07.037, table 1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix) function crystal_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pREAL), dimension(:), intent(in) :: nonSchmidCoefficients !< non-Schmid coefficients for projections real(pREAL), dimension(:), intent(in) :: nonSchmidCoefficients !< non-Schmid coefficients for projections
@ -608,11 +608,11 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
integer :: i integer :: i
if (abs(sense) /= 1) error stop 'Sense in lattice_nonSchmidMatrix' if (abs(sense) /= 1) error stop 'Sense in crystal_nonSchmidMatrix'
coordinateSystem = buildCoordinateSystem(Nslip,CI_NSLIPSYSTEM,CI_SYSTEMSLIP,'cI',0.0_pREAL) coordinateSystem = buildCoordinateSystem(Nslip,CI_NSLIPSYSTEM,CI_SYSTEMSLIP,'cI',0.0_pREAL)
coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip))*real(sense,pREAL) ! convert unidirectional coordinate system coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip))*real(sense,pREAL) ! convert unidirectional coordinate system
nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'cI',0.0_pREAL) ! Schmid contribution nonSchmidMatrix = crystal_SchmidMatrix_slip(Nslip,'cI',0.0_pREAL) ! Schmid contribution
do i = 1,sum(Nslip) do i = 1,sum(Nslip)
direction = coordinateSystem(1:3,1,i) direction = coordinateSystem(1:3,1,i)
@ -635,7 +635,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
+ nonSchmidCoefficients(6) * math_outer(direction, direction) + nonSchmidCoefficients(6) * math_outer(direction, direction)
end do end do
end function lattice_nonSchmidMatrix end function crystal_nonSchmidMatrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -644,7 +644,7 @@ end function lattice_nonSchmidMatrix
!> @details https://doi.org/10.1016/j.actamat.2016.12.040 (cF: Tab S4-1, cI: Tab S5-1) !> @details https://doi.org/10.1016/j.actamat.2016.12.040 (cF: Tab S4-1, cI: Tab S5-1)
!> @details https://doi.org/10.1016/j.ijplas.2014.06.010 (hP: Tab 3b) !> @details https://doi.org/10.1016/j.ijplas.2014.06.010 (hP: Tab 3b)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pREAL), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction real(pREAL), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction
@ -950,19 +950,19 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
interactionTypes = TI_INTERACTIONSLIPSLIP interactionTypes = TI_INTERACTIONSLIPSLIP
NslipMax = TI_NSLIPSYSTEM NslipMax = TI_NSLIPSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_SlipBySlip: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes)
end function lattice_interaction_SlipBySlip end function crystal_interaction_SlipBySlip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Twin-twin interaction matrix !> @brief Twin-twin interaction matrix
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
real(pREAL), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction real(pREAL), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction
@ -1049,19 +1049,19 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(
interactionTypes = HP_INTERACTIONTWINTWIN interactionTypes = HP_INTERACTIONTWINTWIN
NtwinMax = HP_NTWINSYSTEM NtwinMax = HP_NTWINSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_TwinByTwin: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes)
end function lattice_interaction_TwinByTwin end function crystal_interaction_TwinByTwin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Trans-trans interaction matrix !> @brief Trans-trans interaction matrix
!> details only active trans systems are considered !> details only active trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family
real(pREAL), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction real(pREAL), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction
@ -1091,19 +1091,19 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) resu
interactionTypes = CF_INTERACTIONTRANSTRANS interactionTypes = CF_INTERACTIONTRANSTRANS
NtransMax = CF_NTRANSSYSTEM NtransMax = CF_NTRANSSYSTEM
else else
call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_TransByTrans: '//trim(lattice))
end if end if
interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes)
end function lattice_interaction_TransByTrans end function crystal_interaction_TransByTrans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Slip-twin interaction matrix !> @brief Slip-twin interaction matrix
!> details only active slip and twin systems are considered !> details only active slip and twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntwin !< number of active twin systems per family Ntwin !< number of active twin systems per family
@ -1251,19 +1251,19 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
NslipMax = HP_NSLIPSYSTEM NslipMax = HP_NSLIPSYSTEM
NtwinMax = HP_NTWINSYSTEM NtwinMax = HP_NTWINSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_SlipByTwin: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes)
end function lattice_interaction_SlipByTwin end function crystal_interaction_SlipByTwin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Slip-trans interaction matrix !> @brief Slip-trans interaction matrix
!> details only active slip and trans systems are considered !> details only active slip and trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntrans !< number of active trans systems per family Ntrans !< number of active trans systems per family
@ -1304,19 +1304,19 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice)
NslipMax = CF_NSLIPSYSTEM NslipMax = CF_NSLIPSYSTEM
NtransMax = CF_NTRANSSYSTEM NtransMax = CF_NTRANSSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_SlipByTrans: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes)
end function lattice_interaction_SlipByTrans end function crystal_interaction_SlipByTrans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Twin-slip interaction matrix !> @brief Twin-slip interaction matrix
!> details only active twin and slip systems are considered !> details only active twin and slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix) function crystal_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix)
integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family
Nslip !< number of active slip systems per family Nslip !< number of active slip systems per family
@ -1380,19 +1380,19 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) r
NtwinMax = HP_NTWINSYSTEM NtwinMax = HP_NTWINSYSTEM
NslipMax = HP_NSLIPSYSTEM NslipMax = HP_NSLIPSYSTEM
case default case default
call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(lattice)) call IO_error(137,ext_msg='crystal_interaction_TwinBySlip: '//trim(lattice))
end select end select
interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes)
end function lattice_interaction_TwinBySlip end function crystal_interaction_TwinBySlip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Schmid matrix for slip !> @brief Schmid matrix for slip
!> details only active slip systems are considered !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix) function crystal_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1419,7 +1419,7 @@ function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
slipSystems = TI_SYSTEMSLIP slipSystems = TI_SYSTEMSLIP
case default case default
allocate(NslipMax(0)) allocate(NslipMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(lattice)) call IO_error(137,ext_msg='crystal_SchmidMatrix_slip: '//trim(lattice))
end select end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
@ -1435,14 +1435,14 @@ function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
error stop 'dilatational Schmid matrix for slip' error stop 'dilatational Schmid matrix for slip'
end do end do
end function lattice_SchmidMatrix_slip end function crystal_SchmidMatrix_slip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Schmid matrix for twinning !> @brief Schmid matrix for twinning
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix) function crystal_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1466,7 +1466,7 @@ function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
twinSystems = HP_SYSTEMTWIN twinSystems = HP_SYSTEMTWIN
case default case default
allocate(NtwinMax(0)) allocate(NtwinMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_SchmidMatrix_twin: '//trim(lattice))
end select end select
if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
@ -1482,43 +1482,43 @@ function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
error stop 'dilatational Schmid matrix for twin' error stop 'dilatational Schmid matrix for twin'
end do end do
end function lattice_SchmidMatrix_twin end function crystal_SchmidMatrix_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Schmid matrix for transformation !> @brief Schmid matrix for transformation
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_cF,a_cI) result(SchmidMatrix) function crystal_SchmidMatrix_trans(Ntrans,crystal_target,cOverA,a_cF,a_cI) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=*), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: crystal_target !< Bravais lattice (Pearson symbol)
real(pREAL), optional, intent(in) :: cOverA, a_cI, a_cF real(pREAL), optional, intent(in) :: cOverA, a_cI, a_cF
real(pREAL), dimension(3,3,sum(Ntrans)) :: SchmidMatrix real(pREAL), dimension(3,3,sum(Ntrans)) :: SchmidMatrix
real(pREAL), dimension(3,3,sum(Ntrans)) :: devNull real(pREAL), dimension(3,3,sum(Ntrans)) :: devNull
if (lattice_target == 'hP' .and. present(cOverA)) then if (crystal_target == 'hP' .and. present(cOverA)) then
if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) & if (cOverA < 1.0_pREAL .or. cOverA > 2.0_pREAL) &
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) call IO_error(131,ext_msg='crystal_SchmidMatrix_trans: '//trim(crystal_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA=cOverA) call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA=cOverA)
else if (lattice_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then else if (crystal_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then
if (a_cI <= 0.0_pREAL .or. a_cF <= 0.0_pREAL) & if (a_cI <= 0.0_pREAL .or. a_cF <= 0.0_pREAL) &
call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) call IO_error(134,ext_msg='crystal_SchmidMatrix_trans: '//trim(crystal_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,a_cF=a_cF,a_cI=a_cI) call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,a_cF=a_cF,a_cI=a_cI)
else else
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) call IO_error(131,ext_msg='crystal_SchmidMatrix_trans: '//trim(crystal_target))
end if end if
end function lattice_SchmidMatrix_trans end function crystal_SchmidMatrix_trans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Schmid matrix for cleavage !> @brief Schmid matrix for cleavage
!> details only active cleavage systems are considered !> details only active cleavage systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMatrix) function crystal_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1539,7 +1539,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMa
cleavageSystems = CI_SYSTEMCLEAVAGE cleavageSystems = CI_SYSTEMCLEAVAGE
case default case default
allocate(NcleavageMax(0)) allocate(NcleavageMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(lattice)) call IO_error(137,ext_msg='crystal_SchmidMatrix_cleavage: '//trim(lattice))
end select end select
if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0)) & if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0)) &
@ -1555,13 +1555,13 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMa
SchmidMatrix(1:3,1:3,3,i) = math_outer(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,3,i) = math_outer(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,i))
end do end do
end function lattice_SchmidMatrix_cleavage end function crystal_SchmidMatrix_cleavage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Slip direction of slip systems (|| b) !> @brief Slip direction of slip systems (|| b)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_slip_direction(Nslip,lattice,cOverA) result(d) function crystal_slip_direction(Nslip,lattice,cOverA) result(d)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1573,13 +1573,13 @@ function lattice_slip_direction(Nslip,lattice,cOverA) result(d)
coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA) coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
d = coordinateSystem(1:3,1,1:sum(Nslip)) d = coordinateSystem(1:3,1,1:sum(Nslip))
end function lattice_slip_direction end function crystal_slip_direction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Normal direction of slip systems (|| n) !> @brief Normal direction of slip systems (|| n)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_slip_normal(Nslip,lattice,cOverA) result(n) function crystal_slip_normal(Nslip,lattice,cOverA) result(n)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1591,13 +1591,13 @@ function lattice_slip_normal(Nslip,lattice,cOverA) result(n)
coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA) coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
n = coordinateSystem(1:3,2,1:sum(Nslip)) n = coordinateSystem(1:3,2,1:sum(Nslip))
end function lattice_slip_normal end function crystal_slip_normal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Transverse direction of slip systems (|| t = b x n) !> @brief Transverse direction of slip systems (|| t = b x n)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_slip_transverse(Nslip,lattice,cOverA) result(t) function crystal_slip_transverse(Nslip,lattice,cOverA) result(t)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1609,14 +1609,14 @@ function lattice_slip_transverse(Nslip,lattice,cOverA) result(t)
coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA) coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
t = coordinateSystem(1:3,3,1:sum(Nslip)) t = coordinateSystem(1:3,3,1:sum(Nslip))
end function lattice_slip_transverse end function crystal_slip_transverse
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Labels of slip systems !> @brief Labels of slip systems
!> details only active slip systems are considered !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_labels_slip(Nslip,lattice) result(labels) function crystal_labels_slip(Nslip,lattice) result(labels)
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1640,7 +1640,7 @@ function lattice_labels_slip(Nslip,lattice) result(labels)
NslipMax = TI_NSLIPSYSTEM NslipMax = TI_NSLIPSYSTEM
slipSystems = TI_SYSTEMSLIP slipSystems = TI_SYSTEMSLIP
case default case default
call IO_error(137,ext_msg='lattice_labels_slip: '//trim(lattice)) call IO_error(137,ext_msg='crystal_labels_slip: '//trim(lattice))
end select end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
@ -1650,13 +1650,13 @@ function lattice_labels_slip(Nslip,lattice) result(labels)
labels = getLabels(Nslip,NslipMax,slipSystems) labels = getLabels(Nslip,NslipMax,slipSystems)
end function lattice_labels_slip end function crystal_labels_slip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return 3x3 tensor with symmetry according to given Bravais lattice !> @brief Return 3x3 tensor with symmetry according to given Bravais lattice
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_symmetrize_33(T,lattice) result(T_sym) pure function crystal_symmetrize_33(T,lattice) result(T_sym)
real(pREAL), dimension(3,3) :: T_sym real(pREAL), dimension(3,3) :: T_sym
@ -1677,14 +1677,14 @@ pure function lattice_symmetrize_33(T,lattice) result(T_sym)
T_sym(3,3) = T(3,3) T_sym(3,3) = T(3,3)
end select end select
end function lattice_symmetrize_33 end function crystal_symmetrize_33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice !> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice
!> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym) pure function crystal_symmetrize_C66(C66,lattice) result(C66_sym)
real(pREAL), dimension(6,6) :: C66_sym real(pREAL), dimension(6,6) :: C66_sym
@ -1723,14 +1723,14 @@ pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym)
end do end do
end do end do
end function lattice_symmetrize_C66 end function crystal_symmetrize_C66
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Labels for twin systems !> @brief Labels for twin systems
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_labels_twin(Ntwin,lattice) result(labels) function crystal_labels_twin(Ntwin,lattice) result(labels)
integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family
character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol)
@ -1751,7 +1751,7 @@ function lattice_labels_twin(Ntwin,lattice) result(labels)
NtwinMax = HP_NTWINSYSTEM NtwinMax = HP_NTWINSYSTEM
twinSystems = HP_SYSTEMTWIN twinSystems = HP_SYSTEMTWIN
case default case default
call IO_error(137,ext_msg='lattice_labels_twin: '//trim(lattice)) call IO_error(137,ext_msg='crystal_labels_twin: '//trim(lattice))
end select end select
if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
@ -1761,7 +1761,7 @@ function lattice_labels_twin(Ntwin,lattice) result(labels)
labels = getLabels(Ntwin,NtwinMax,twinSystems) labels = getLabels(Ntwin,NtwinMax,twinSystems)
end function lattice_labels_twin end function crystal_labels_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1778,8 +1778,8 @@ function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
real(pREAL), dimension(3,sum(Nslip)) :: n, t real(pREAL), dimension(3,sum(Nslip)) :: n, t
integer :: i, j integer :: i, j
n = lattice_slip_normal (Nslip,lattice,cOverA) n = crystal_slip_normal (Nslip,lattice,cOverA)
t = lattice_slip_transverse(Nslip,lattice,cOverA) t = crystal_slip_transverse(Nslip,lattice,cOverA)
do i=1, sum(Nslip); do j=1, sum(Nslip) do i=1, sum(Nslip); do j=1, sum(Nslip)
projection(i,j) = abs(math_inner(n(:,i),t(:,j))) projection(i,j) = abs(math_inner(n(:,i),t(:,j)))
@ -1802,8 +1802,8 @@ function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
real(pREAL), dimension(3,sum(Nslip)) :: n, d real(pREAL), dimension(3,sum(Nslip)) :: n, d
integer :: i, j integer :: i, j
n = lattice_slip_normal (Nslip,lattice,cOverA) n = crystal_slip_normal (Nslip,lattice,cOverA)
d = lattice_slip_direction(Nslip,lattice,cOverA) d = crystal_slip_direction(Nslip,lattice,cOverA)
do i=1, sum(Nslip); do j=1, sum(Nslip) do i=1, sum(Nslip); do j=1, sum(Nslip)
projection(i,j) = abs(math_inner(n(:,i),d(:,j))) projection(i,j) = abs(math_inner(n(:,i),d(:,j)))
@ -2150,7 +2150,7 @@ end function getlabels
!> @brief Equivalent Poisson's ratio (ν) !> @brief Equivalent Poisson's ratio (ν)
!> @details https://doi.org/10.1143/JPSJ.20.635 !> @details https://doi.org/10.1143/JPSJ.20.635
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_isotropic_nu(C,assumption,lattice) result(nu) pure function crystal_isotropic_nu(C,assumption,lattice) result(nu)
real(pREAL), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pREAL), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss') character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss')
@ -2172,10 +2172,10 @@ pure function lattice_isotropic_nu(C,assumption,lattice) result(nu)
error stop 'invalid assumption' error stop 'invalid assumption'
end if end if
mu = lattice_isotropic_mu(C,assumption,lattice) mu = crystal_isotropic_mu(C,assumption,lattice)
nu = (1.5_pREAL*K-mu)/(3.0_pREAL*K+mu) nu = (1.5_pREAL*K-mu)/(3.0_pREAL*K+mu)
end function lattice_isotropic_nu end function crystal_isotropic_nu
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2183,7 +2183,7 @@ end function lattice_isotropic_nu
!> @details https://doi.org/10.1143/JPSJ.20.635 !> @details https://doi.org/10.1143/JPSJ.20.635
!> @details Nonlinear Mechanics of Crystals 10.1007/978-94-007-0350-6, pp 563 !> @details Nonlinear Mechanics of Crystals 10.1007/978-94-007-0350-6, pp 563
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_isotropic_mu(C,assumption,lattice) result(mu) pure function crystal_isotropic_mu(C,assumption,lattice) result(mu)
real(pREAL), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) real(pREAL), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss') character(len=*), intent(in) :: assumption !< Assumption (isostrain = 'Voigt', isostress = 'Reuss')
@ -2220,11 +2220,11 @@ pure function lattice_isotropic_mu(C,assumption,lattice) result(mu)
error stop 'invalid assumption' error stop 'invalid assumption'
end if end if
end function lattice_isotropic_mu end function crystal_isotropic_mu
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Check correctness of some lattice functions. !> @brief Check correctness of some crystal functions.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine selfTest subroutine selfTest
@ -2246,10 +2246,10 @@ subroutine selfTest
do i = 1, 10 do i = 1, 10
call random_number(C) call random_number(C)
C_cF = lattice_symmetrize_C66(C,'cI') C_cF = crystal_symmetrize_C66(C,'cI')
C_cI = lattice_symmetrize_C66(C,'cF') C_cI = crystal_symmetrize_C66(C,'cF')
C_hP = lattice_symmetrize_C66(C,'hP') C_hP = crystal_symmetrize_C66(C,'hP')
C_tI = lattice_symmetrize_C66(C,'tI') C_tI = crystal_symmetrize_C66(C,'tI')
if (any(dNeq(C_cI,transpose(C_cF)))) error stop 'SymmetryC66/cI-cF' if (any(dNeq(C_cI,transpose(C_cF)))) error stop 'SymmetryC66/cI-cF'
if (any(dNeq(C_cF,transpose(C_cI)))) error stop 'SymmetryC66/cF-cI' if (any(dNeq(C_cF,transpose(C_cI)))) error stop 'SymmetryC66/cF-cI'
@ -2269,10 +2269,10 @@ subroutine selfTest
if (any(dNeq(C(4,4),[C_tI(4,4),C_tI(5,5)]))) error stop 'SymmetryC_44-55/tI' if (any(dNeq(C(4,4),[C_tI(4,4),C_tI(5,5)]))) error stop 'SymmetryC_44-55/tI'
call random_number(T) call random_number(T)
T_cF = lattice_symmetrize_33(T,'cI') T_cF = crystal_symmetrize_33(T,'cI')
T_cI = lattice_symmetrize_33(T,'cF') T_cI = crystal_symmetrize_33(T,'cF')
T_hP = lattice_symmetrize_33(T,'hP') T_hP = crystal_symmetrize_33(T,'hP')
T_tI = lattice_symmetrize_33(T,'tI') T_tI = crystal_symmetrize_33(T,'tI')
if (any(dNeq0(T_cF) .and. math_I3<1.0_pREAL)) error stop 'Symmetry33/c' if (any(dNeq0(T_cF) .and. math_I3<1.0_pREAL)) error stop 'Symmetry33/c'
if (any(dNeq0(T_hP) .and. math_I3<1.0_pREAL)) error stop 'Symmetry33/hP' if (any(dNeq0(T_hP) .and. math_I3<1.0_pREAL)) error stop 'Symmetry33/hP'
@ -2291,48 +2291,48 @@ subroutine selfTest
C(4,4) = 0.5_pREAL * (C(1,1) - C(1,2)) C(4,4) = 0.5_pREAL * (C(1,1) - C(1,2))
C(6,6) = C(4,4) C(6,6) = C(4,4)
C_cI = lattice_symmetrize_C66(C,'cI') C_cI = crystal_symmetrize_C66(C,'cI')
if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'isostrain','cI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/cI' if (dNeq(C_cI(4,4),crystal_isotropic_mu(C_cI,'isostrain','cI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/cI'
if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'isostress','cI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/cI' if (dNeq(C_cI(4,4),crystal_isotropic_mu(C_cI,'isostress','cI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/cI'
lambda = C_cI(1,2) lambda = C_cI(1,2)
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_cI,'isostrain','cI')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_cI,'isostrain','cI')), &
lattice_isotropic_nu(C_cI,'isostrain','cI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/cI' crystal_isotropic_nu(C_cI,'isostrain','cI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/cI'
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_cI,'isostress','cI')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_cI,'isostress','cI')), &
lattice_isotropic_nu(C_cI,'isostress','cI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/cI' crystal_isotropic_nu(C_cI,'isostress','cI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/cI'
C_hP = lattice_symmetrize_C66(C,'hP') C_hP = crystal_symmetrize_C66(C,'hP')
if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'isostrain','hP'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/hP' if (dNeq(C(4,4),crystal_isotropic_mu(C_hP,'isostrain','hP'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/hP'
if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'isostress','hP'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/hP' if (dNeq(C(4,4),crystal_isotropic_mu(C_hP,'isostress','hP'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/hP'
lambda = C_hP(1,2) lambda = C_hP(1,2)
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_hP,'isostrain','hP')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_hP,'isostrain','hP')), &
lattice_isotropic_nu(C_hP,'isostrain','hP'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/hP' crystal_isotropic_nu(C_hP,'isostrain','hP'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/hP'
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_hP,'isostress','hP')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_hP,'isostress','hP')), &
lattice_isotropic_nu(C_hP,'isostress','hP'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/hP' crystal_isotropic_nu(C_hP,'isostress','hP'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/hP'
C_tI = lattice_symmetrize_C66(C,'tI') C_tI = crystal_symmetrize_C66(C,'tI')
if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'isostrain','tI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/tI' if (dNeq(C(6,6),crystal_isotropic_mu(C_tI,'isostrain','tI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostrain/tI'
if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'isostress','tI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/tI' if (dNeq(C(6,6),crystal_isotropic_mu(C_tI,'isostress','tI'),1.0e-12_pREAL)) error stop 'isotropic_mu/isostress/tI'
lambda = C_tI(1,2) lambda = C_tI(1,2)
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_tI,'isostrain','tI')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_tI,'isostrain','tI')), &
lattice_isotropic_nu(C_tI,'isostrain','tI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/tI' crystal_isotropic_nu(C_tI,'isostrain','tI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostrain/tI'
if (dNeq(lambda*0.5_pREAL/(lambda+lattice_isotropic_mu(C_tI,'isostress','tI')), & if (dNeq(lambda*0.5_pREAL/(lambda+crystal_isotropic_mu(C_tI,'isostress','tI')), &
lattice_isotropic_nu(C_tI,'isostress','tI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/tI' crystal_isotropic_nu(C_tI,'isostress','tI'),1.0e-12_pREAL)) error stop 'isotropic_nu/isostress/tI'
call random_number(C) call random_number(C)
C = lattice_symmetrize_C66(C+math_eye(6),'cI') C = crystal_symmetrize_C66(C+math_eye(6),'cI')
if (dNeq(lattice_isotropic_mu(C,'isostrain','cI'), lattice_isotropic_mu(C,'isostrain','hP'), 1.0e-12_pREAL)) & if (dNeq(crystal_isotropic_mu(C,'isostrain','cI'), crystal_isotropic_mu(C,'isostrain','hP'), 1.0e-12_pREAL)) &
error stop 'isotropic_mu/isostrain/cI-hP' error stop 'isotropic_mu/isostrain/cI-hP'
if (dNeq(lattice_isotropic_nu(C,'isostrain','cF'), lattice_isotropic_nu(C,'isostrain','cI'), 1.0e-12_pREAL)) & if (dNeq(crystal_isotropic_nu(C,'isostrain','cF'), crystal_isotropic_nu(C,'isostrain','cI'), 1.0e-12_pREAL)) &
error stop 'isotropic_nu/isostrain/cF-tI' error stop 'isotropic_nu/isostrain/cF-tI'
if (dNeq(lattice_isotropic_mu(C,'isostress','cI'), lattice_isotropic_mu(C,'isostress'), 1.0e-12_pREAL)) & if (dNeq(crystal_isotropic_mu(C,'isostress','cI'), crystal_isotropic_mu(C,'isostress'), 1.0e-12_pREAL)) &
error stop 'isotropic_mu/isostress/cI-hP' error stop 'isotropic_mu/isostress/cI-hP'
if (dNeq(lattice_isotropic_nu(C,'isostress','cF'), lattice_isotropic_nu(C,'isostress'), 1.0e-12_pREAL)) & if (dNeq(crystal_isotropic_nu(C,'isostress','cF'), crystal_isotropic_nu(C,'isostress'), 1.0e-12_pREAL)) &
error stop 'isotropic_nu/isostress/cF-tI' error stop 'isotropic_nu/isostress/cF-tI'
end subroutine selfTest end subroutine selfTest
end module lattice end module crystal

44
src/grid/DAMASK_grid.f90
View File

@ -24,7 +24,7 @@ program DAMASK_grid
use material use material
use spectral_utilities use spectral_utilities
use grid_mechanical_spectral_basic use grid_mechanical_spectral_basic
use grid_mechanical_spectral_polarisation use grid_mechanical_spectral_polarization
use grid_mechanical_FEM use grid_mechanical_FEM
use grid_damage_spectral use grid_damage_spectral
use grid_thermal_spectral use grid_thermal_spectral
@ -107,7 +107,8 @@ program DAMASK_grid
external :: & external :: &
quit quit
type(tDict), pointer :: & type(tDict), pointer :: &
config_load, & load, &
num_solver, &
num_grid, & num_grid, &
load_step, & load_step, &
solver, & solver, &
@ -122,6 +123,7 @@ program DAMASK_grid
character(len=:), allocatable :: & character(len=:), allocatable :: &
fileContent, fname fileContent, fname
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
@ -134,12 +136,16 @@ program DAMASK_grid
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read (and check) field parameters from numerics file ! read (and check) field parameters from numerics file
num_grid => config_numerics%get_dict('grid', defaultVal=emptyDict)
stagItMax = num_grid%get_asInt('maxStaggeredIter',defaultVal=10)
maxCutBack = num_grid%get_asInt('maxCutBack',defaultVal=3)
if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict)
if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') num_grid => num_solver%get_dict('grid',defaultVal=emptyDict)
stagItMax = num_grid%get_asInt('N_staggered_iter_max',defaultVal=10)
maxCutBack = num_grid%get_asInt('N_cutback_max',defaultVal=3)
if (stagItMax < 0) call IO_error(301,ext_msg='N_staggered_iter_max')
if (maxCutBack < 0) call IO_error(301,ext_msg='N_cutback_max')
if (worldrank == 0) then if (worldrank == 0) then
fileContent = IO_read(CLI_loadFile) fileContent = IO_read(CLI_loadFile)
@ -151,8 +157,8 @@ program DAMASK_grid
end if end if
call parallelization_bcast_str(fileContent) call parallelization_bcast_str(fileContent)
config_load => YAML_parse_str_asDict(fileContent) load => YAML_parse_str_asDict(fileContent)
solver => config_load%get_dict('solver') solver => load%get_dict('solver')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assign mechanics solver depending on selected type ! assign mechanics solver depending on selected type
@ -167,11 +173,11 @@ program DAMASK_grid
mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite
case ('spectral_polarization') case ('spectral_polarization')
mechanical_init => grid_mechanical_spectral_polarisation_init mechanical_init => grid_mechanical_spectral_polarization_init
mechanical_forward => grid_mechanical_spectral_polarisation_forward mechanical_forward => grid_mechanical_spectral_polarization_forward
mechanical_solution => grid_mechanical_spectral_polarisation_solution mechanical_solution => grid_mechanical_spectral_polarization_solution
mechanical_updateCoords => grid_mechanical_spectral_polarisation_updateCoords mechanical_updateCoords => grid_mechanical_spectral_polarization_updateCoords
mechanical_restartWrite => grid_mechanical_spectral_polarisation_restartWrite mechanical_restartWrite => grid_mechanical_spectral_polarization_restartWrite
case ('FEM') case ('FEM')
mechanical_init => grid_mechanical_FEM_init mechanical_init => grid_mechanical_FEM_init
@ -206,7 +212,7 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
load_steps => config_load%get_list('loadstep') load_steps => load%get_list('loadstep')
allocate(loadCases(load_steps%length)) ! array of load cases allocate(loadCases(load_steps%length)) ! array of load cases
do l = 1, load_steps%length do l = 1, load_steps%length
@ -313,21 +319,21 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers ! doing initialization depending on active solvers
call spectral_Utilities_init() call spectral_utilities_init()
do field = 2, nActiveFields do field = 2, nActiveFields
select case (ID(field)) select case (ID(field))
case (FIELD_THERMAL_ID) case (FIELD_THERMAL_ID)
call grid_thermal_spectral_init() call grid_thermal_spectral_init(num_grid)
case (FIELD_DAMAGE_ID) case (FIELD_DAMAGE_ID)
call grid_damage_spectral_init() call grid_damage_spectral_init(num_grid)
end select end select
end do end do
call mechanical_init() call mechanical_init(num_grid)
call config_numerics_deallocate() call config_numerics_deallocate()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

43
src/grid/grid_damage_spectral.f90
View File

@ -16,6 +16,7 @@ module grid_damage_spectral
use prec use prec
use parallelization use parallelization
use IO use IO
use misc
use CLI use CLI
use HDF5_utilities use HDF5_utilities
use HDF5 use HDF5
@ -67,9 +68,11 @@ module grid_damage_spectral
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data !> @brief Allocate all necessary fields and fill them with data, potentially from restart file.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_init() subroutine grid_damage_spectral_init(num_grid)
type(tDict), pointer, intent(in) :: num_grid
integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global
integer :: i, j, k, ce integer :: i, j, k, ce
@ -81,10 +84,13 @@ subroutine grid_damage_spectral_init()
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
real(pREAL), dimension(1,product(cells(1:2))*cells3) :: tempN real(pREAL), dimension(1,product(cells(1:2))*cells3) :: tempN
type(tDict), pointer :: & type(tDict), pointer :: &
num_grid, & num_grid_damage
num_generic
character(len=pSTRLEN) :: & character(len=pSTRLEN) :: &
snes_type snes_type
character(len=:), allocatable :: &
extmsg, &
petsc_options
print'(/,1x,a)', '<<<+- grid_spectral_damage init -+>>>' print'(/,1x,a)', '<<<+- grid_spectral_damage init -+>>>'
@ -96,25 +102,26 @@ subroutine grid_damage_spectral_init()
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict) num_grid_damage => num_grid%get_dict('damage',defaultVal=emptyDict)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
num%eps_damage_atol = num_grid%get_asReal ('eps_damage_atol',defaultVal=1.0e-2_pREAL)
num%eps_damage_rtol = num_grid%get_asReal ('eps_damage_rtol',defaultVal=1.0e-6_pREAL)
num_generic => config_numerics%get_dict('generic',defaultVal=emptyDict) num%itmax = num_grid_damage%get_asInt('N_iter_max', defaultVal=100)
num%phi_min = num_generic%get_asReal('phi_min', defaultVal=1.0e-6_pREAL) num%eps_damage_atol = num_grid_damage%get_asReal('eps_abs_phi',defaultVal=1.0e-2_pREAL)
num%eps_damage_rtol = num_grid_damage%get_asReal('eps_rel_phi',defaultVal=1.0e-6_pREAL)
num%phi_min = num_grid_damage%get_asReal('phi_min', defaultVal=1.0e-6_pREAL)
if (num%phi_min < 0.0_pREAL) call IO_error(301,ext_msg='phi_min') extmsg = ''
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%eps_damage_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_phi'
if (num%eps_damage_atol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_damage_atol') if (num%eps_damage_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_phi'
if (num%eps_damage_rtol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_damage_rtol') if (num%phi_min < 0.0_pREAL) extmsg = trim(extmsg)//' phi_min'
if (num%itmax < 1) extmsg = trim(extmsg)//' N_iter_max'
if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-damage_snes_type newtonls -damage_snes_mf & petsc_options = misc_prefixOptions('-snes_type newtonls -snes_mf -snes_ksp_ew -ksp_type fgmres '// &
&-damage_snes_ksp_ew -damage_ksp_type fgmres',err_PETSc) num_grid_damage%get_asStr('PETSc_options',defaultVal=''),'damage_')
CHKERRQ(err_PETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asStr('petsc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

53
src/grid/grid_mech_FEM.f90
View File

@ -15,8 +15,9 @@ module grid_mechanical_FEM
use prec use prec
use parallelization use parallelization
use CLI
use IO use IO
use misc
use CLI
use HDF5 use HDF5
use HDF5_utilities use HDF5_utilities
use math use math
@ -94,9 +95,11 @@ module grid_mechanical_FEM
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief Allocate all necessary fields and fill them with data, potentially from restart info.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_FEM_init subroutine grid_mechanical_FEM_init(num_grid)
type(tDict), pointer, intent(in) :: num_grid
real(pREAL), parameter :: HGCoeff = 0.0e-2_pREAL real(pREAL), parameter :: HGCoeff = 0.0e-2_pREAL
real(pREAL), parameter, dimension(4,8) :: & real(pREAL), parameter, dimension(4,8) :: &
@ -118,41 +121,41 @@ subroutine grid_mechanical_FEM_init
integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
type(tDict), pointer :: & type(tDict), pointer :: &
num_grid num_grid_mech
character(len=pSTRLEN) :: & character(len=:), allocatable :: &
extmsg = '' extmsg, &
petsc_options
print'(/,1x,a)', '<<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT)
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict) num_grid_mech => num_grid%get_dict('mechanical',defaultVal=emptyDict)
num%eps_div_atol = num_grid%get_asReal('eps_div_atol', defaultVal=1.0e-4_pREAL) num%itmin = num_grid_mech%get_asInt('N_iter_min',defaultVal=1)
num%eps_div_rtol = num_grid%get_asReal('eps_div_rtol', defaultVal=5.0e-4_pREAL) num%itmax = num_grid_mech%get_asInt('N_iter_max',defaultVal=100)
num%eps_stress_atol = num_grid%get_asReal('eps_stress_atol',defaultVal=1.0e3_pREAL) num%eps_div_atol = num_grid_mech%get_asReal('eps_abs_div(P)',defaultVal=1.0e-4_pREAL)
num%eps_stress_rtol = num_grid%get_asReal('eps_stress_rtol',defaultVal=1.0e-3_pREAL) num%eps_div_rtol = num_grid_mech%get_asReal('eps_rel_div(P)',defaultVal=5.0e-4_pREAL)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) num%eps_stress_atol = num_grid_mech%get_asReal('eps_abs_P', defaultVal=1.0e3_pREAL)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) num%eps_stress_rtol = num_grid_mech%get_asReal('eps_rel_P', defaultVal=1.0e-3_pREAL)
if (num%eps_div_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_div_atol' extmsg = ''
if (num%eps_div_rtol < 0.0_pREAL) extmsg = trim(extmsg)//' eps_div_rtol' if (num%eps_div_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_div(P)'
if (num%eps_stress_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_stress_atol' if (num%eps_div_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_div(P)'
if (num%eps_stress_rtol < 0.0_pREAL) extmsg = trim(extmsg)//' eps_stress_rtol' if (num%eps_stress_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_P'
if (num%itmax <= 1) extmsg = trim(extmsg)//' itmax' if (num%eps_stress_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_P'
if (num%itmin > num%itmax .or. num%itmin < 1) extmsg = trim(extmsg)//' itmin' if (num%itmax < 1) extmsg = trim(extmsg)//' N_iter_max'
if (num%itmin > num%itmax .or. num%itmin < 1) extmsg = trim(extmsg)//' N_iter_min'
if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg)) if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS, &
'-mechanical_snes_type newtonls -mechanical_ksp_type fgmres & petsc_options = misc_prefixOptions('-snes_type newtonls -ksp_type fgmres -ksp_max_it 25 '// &
&-mechanical_ksp_max_it 25', & num_grid_mech%get_asStr('PETSc_options',defaultVal='') ,'mechanical_')
err_PETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,err_PETSc)
CHKERRQ(err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asStr('petsc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

52
src/grid/grid_mech_spectral_basic.f90
View File

@ -16,13 +16,13 @@ module grid_mechanical_spectral_basic
use prec use prec
use parallelization use parallelization
use CLI use CLI
use misc
use IO use IO
use HDF5 use HDF5
use HDF5_utilities use HDF5_utilities
use math use math
use rotations use rotations
use spectral_utilities use spectral_utilities
use config
use homogenization use homogenization
use discretization_grid use discretization_grid
@ -101,9 +101,11 @@ module grid_mechanical_spectral_basic
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief Allocate all necessary fields and fill them with data, potentially from restart info.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_basic_init() subroutine grid_mechanical_spectral_basic_init(num_grid)
type(tDict), pointer, intent(in) :: num_grid
real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: P real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: P
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
@ -114,9 +116,11 @@ subroutine grid_mechanical_spectral_basic_init()
real(pREAL), dimension(3,3,product(cells(1:2))*cells3) :: temp33n real(pREAL), dimension(3,3,product(cells(1:2))*cells3) :: temp33n
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
type(tDict), pointer :: & type(tDict), pointer :: &
num_grid num_grid_fft, &
character(len=pSTRLEN) :: & num_grid_mech
extmsg = '' character(len=:), allocatable :: &
extmsg, &
petsc_options
print'(/,1x,a)', '<<<+- grid_mechanical_spectral_basic init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- grid_mechanical_spectral_basic init -+>>>'; flush(IO_STDOUT)
@ -129,30 +133,32 @@ subroutine grid_mechanical_spectral_basic_init()
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict) num_grid_fft => num_grid%get_dict('FFT',defaultVal=emptyDict)
num_grid_mech => num_grid%get_dict('mechanical',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool('update_gamma', defaultVal=.false.) num%itmin = num_grid_mech%get_asInt('N_iter_min',defaultVal=1)
num%eps_div_atol = num_grid%get_asReal('eps_div_atol', defaultVal=1.0e-4_pREAL) num%itmax = num_grid_mech%get_asInt('N_iter_max',defaultVal=100)
num%eps_div_rtol = num_grid%get_asReal('eps_div_rtol', defaultVal=5.0e-4_pREAL) num%update_gamma = num_grid_mech%get_asBool('update_gamma',defaultVal=.false.)
num%eps_stress_atol = num_grid%get_asReal('eps_stress_atol',defaultVal=1.0e3_pREAL) num%eps_div_atol = num_grid_mech%get_asReal('eps_abs_div(P)', defaultVal=1.0e-4_pREAL)
num%eps_stress_rtol = num_grid%get_asReal('eps_stress_rtol',defaultVal=1.0e-3_pREAL) num%eps_div_rtol = num_grid_mech%get_asReal('eps_rel_div(P)', defaultVal=5.0e-4_pREAL)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) num%eps_stress_atol = num_grid_mech%get_asReal('eps_abs_P', defaultVal=1.0e3_pREAL)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) num%eps_stress_rtol = num_grid_mech%get_asReal('eps_rel_P', defaultVal=1.0e-3_pREAL)
if (num%eps_div_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_div_atol' extmsg = ''
if (num%eps_div_rtol < 0.0_pREAL) extmsg = trim(extmsg)//' eps_div_rtol' if (num%eps_div_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_div(P)'
if (num%eps_stress_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_stress_atol' if (num%eps_div_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_div(P)'
if (num%eps_stress_rtol < 0.0_pREAL) extmsg = trim(extmsg)//' eps_stress_rtol' if (num%eps_stress_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_P'
if (num%itmax <= 1) extmsg = trim(extmsg)//' itmax' if (num%eps_stress_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_P'
if (num%itmin > num%itmax .or. num%itmin < 1) extmsg = trim(extmsg)//' itmin' if (num%itmax < 1) extmsg = trim(extmsg)//' N_iter_max'
if (num%itmin > num%itmax .or. num%itmin < 1) extmsg = trim(extmsg)//' N_iter_min'
if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg)) if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type ngmres',err_PETSc) petsc_options = misc_prefixOptions('-snes_type ngmres '//num_grid_mech%get_asStr('PETSc_options',defaultVal=''), &
CHKERRQ(err_PETSc) 'mechanical_')
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asStr('petsc_options',defaultVal=''),err_PETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

96
src/grid/grid_mech_spectral_polarisation.f90 → src/grid/grid_mech_spectral_polarization.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Grid solver for mechanics: Spectral Polarisation !> @brief Grid solver for mechanics: Spectral Polarisation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module grid_mechanical_spectral_polarisation module grid_mechanical_spectral_polarization
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScDMDA use PETScDMDA
@ -16,6 +16,7 @@ module grid_mechanical_spectral_polarisation
use prec use prec
use parallelization use parallelization
use CLI use CLI
use misc
use IO use IO
use HDF5 use HDF5
use HDF5_utilities use HDF5_utilities
@ -103,18 +104,20 @@ module grid_mechanical_spectral_polarisation
totalIter = 0 !< total iteration in current increment totalIter = 0 !< total iteration in current increment
public :: & public :: &
grid_mechanical_spectral_polarisation_init, & grid_mechanical_spectral_polarization_init, &
grid_mechanical_spectral_polarisation_solution, & grid_mechanical_spectral_polarization_solution, &
grid_mechanical_spectral_polarisation_forward, & grid_mechanical_spectral_polarization_forward, &
grid_mechanical_spectral_polarisation_updateCoords, & grid_mechanical_spectral_polarization_updateCoords, &
grid_mechanical_spectral_polarisation_restartWrite grid_mechanical_spectral_polarization_restartWrite
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief Allocate all necessary fields and fill them with data, potentially from restart info.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_polarisation_init() subroutine grid_mechanical_spectral_polarization_init(num_grid)
type(tDict), pointer, intent(in) :: num_grid
real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: P real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: P
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
@ -127,9 +130,11 @@ subroutine grid_mechanical_spectral_polarisation_init()
real(pREAL), dimension(3,3,product(cells(1:2))*cells3) :: temp33n real(pREAL), dimension(3,3,product(cells(1:2))*cells3) :: temp33n
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
type(tDict), pointer :: & type(tDict), pointer :: &
num_grid num_grid_fft,&
character(len=pSTRLEN) :: & num_grid_mech
extmsg = '' character(len=:), allocatable :: &
extmsg, &
petsc_options
print '(/,1x,a)', '<<<+- grid_mechanical_spectral_polarization init -+>>>'; flush(IO_STDOUT) print '(/,1x,a)', '<<<+- grid_mechanical_spectral_polarization init -+>>>'; flush(IO_STDOUT)
@ -137,41 +142,42 @@ subroutine grid_mechanical_spectral_polarisation_init()
print '(/,1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print '(/,1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print '( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006' print '( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict) num_grid_fft => num_grid%get_dict('FFT',defaultVal=emptyDict)
num_grid_mech => num_grid%get_dict('mechanical',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool('update_gamma', defaultVal=.false.) num%itmin = num_grid_mech%get_asInt('N_iter_min',defaultVal=1)
num%eps_div_atol = num_grid%get_asReal('eps_div_atol', defaultVal=1.0e-4_pREAL) num%itmax = num_grid_mech%get_asInt('N_iter_max',defaultVal=100)
num%eps_div_rtol = num_grid%get_asReal('eps_div_rtol', defaultVal=5.0e-4_pREAL) num%update_gamma = num_grid_mech%get_asBool('update_gamma',defaultVal=.false.)
num%eps_curl_atol = num_grid%get_asReal('eps_curl_atol', defaultVal=1.0e-10_pREAL) num%eps_div_atol = num_grid_mech%get_asReal('eps_abs_div(P)', defaultVal=1.0e-4_pREAL)
num%eps_curl_rtol = num_grid%get_asReal('eps_curl_rtol', defaultVal=5.0e-4_pREAL) num%eps_div_rtol = num_grid_mech%get_asReal('eps_rel_div(P)', defaultVal=5.0e-4_pREAL)
num%eps_stress_atol = num_grid%get_asReal('eps_stress_atol',defaultVal=1.0e3_pREAL) num%eps_curl_atol = num_grid_mech%get_asReal('eps_abs_curl(F)',defaultVal=1.0e-10_pREAL)
num%eps_stress_rtol = num_grid%get_asReal('eps_stress_rtol',defaultVal=1.0e-3_pREAL) num%eps_curl_rtol = num_grid_mech%get_asReal('eps_rel_curl(F)',defaultVal=5.0e-4_pREAL)
num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) num%eps_stress_atol = num_grid_mech%get_asReal('eps_abs_P', defaultVal=1.0e3_pREAL)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) num%eps_stress_rtol = num_grid_mech%get_asReal('eps_rel_P', defaultVal=1.0e-3_pREAL)
num%alpha = num_grid%get_asReal('alpha', defaultVal=1.0_pREAL) num%alpha = num_grid_mech%get_asReal('alpha', defaultVal=1.0_pREAL)
num%beta = num_grid%get_asReal('beta', defaultVal=1.0_pREAL) num%beta = num_grid_mech%get_asReal('beta', defaultVal=1.0_pREAL)
if (num%eps_div_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_div_atol' extmsg = ''
if (num%eps_div_rtol < 0.0_pREAL) extmsg = trim(extmsg)//' eps_div_rtol' if (num%eps_div_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_div(P)'
if (num%eps_curl_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_curl_atol' if (num%eps_div_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_div(P)'
if (num%eps_curl_rtol < 0.0_pREAL) extmsg = trim(extmsg)//' eps_curl_rtol' if (num%eps_curl_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_curl(F)'
if (num%eps_stress_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_stress_atol' if (num%eps_curl_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_curl(F)'
if (num%eps_stress_rtol < 0.0_pREAL) extmsg = trim(extmsg)//' eps_stress_rtol' if (num%eps_stress_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_P'
if (num%itmax <= 1) extmsg = trim(extmsg)//' itmax' if (num%eps_stress_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_P'
if (num%itmin > num%itmax .or. num%itmin < 1) extmsg = trim(extmsg)//' itmin'
if (num%alpha <= 0.0_pREAL .or. num%alpha > 2.0_pREAL) extmsg = trim(extmsg)//' alpha' if (num%alpha <= 0.0_pREAL .or. num%alpha > 2.0_pREAL) extmsg = trim(extmsg)//' alpha'
if (num%beta < 0.0_pREAL .or. num%beta > 2.0_pREAL) extmsg = trim(extmsg)//' beta' if (num%beta < 0.0_pREAL .or. num%beta > 2.0_pREAL) extmsg = trim(extmsg)//' beta'
if (num%itmax < 1) extmsg = trim(extmsg)//' N_iter_max'
if (num%itmin > num%itmax .or. num%itmin < 1) extmsg = trim(extmsg)//' N_iter_min'
if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg)) if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type ngmres',err_PETSc) petsc_options = misc_prefixOptions('-snes_type ngmres '//num_grid_mech%get_asStr('PETSc_options',defaultVal=''), &
CHKERRQ(err_PETSc) 'mechanical_')
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asStr('petsc_options',defaultVal=''),err_PETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -285,13 +291,13 @@ subroutine grid_mechanical_spectral_polarisation_init()
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg) S_scale = math_invSym3333(C_minMaxAvg)
end subroutine grid_mechanical_spectral_polarisation_init end subroutine grid_mechanical_spectral_polarization_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the Polarisation scheme with internal iterations !> @brief solution for the Polarisation scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(solution) function grid_mechanical_spectral_polarization_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
@ -326,14 +332,14 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_av,P_aim,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_spectral_polarisation_solution end function grid_mechanical_spectral_polarization_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep !> @details find new boundary conditions and best F estimate for end of current timestep
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,& subroutine grid_mechanical_spectral_polarization_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC) deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: & logical, intent(in) :: &
@ -427,13 +433,13 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%Delta_t = Delta_t params%Delta_t = Delta_t
end subroutine grid_mechanical_spectral_polarisation_forward end subroutine grid_mechanical_spectral_polarization_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Update coordinates. !> @brief Update coordinates.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_polarisation_updateCoords() subroutine grid_mechanical_spectral_polarization_updateCoords()
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
real(pREAL), dimension(:,:,:,:), pointer :: FandF_tau real(pREAL), dimension(:,:,:,:), pointer :: FandF_tau
@ -444,13 +450,13 @@ subroutine grid_mechanical_spectral_polarisation_updateCoords()
call DMDAVecRestoreArrayReadF90(DM_mech,FandF_tau_PETSc,FandF_tau,err_PETSc) call DMDAVecRestoreArrayReadF90(DM_mech,FandF_tau_PETSc,FandF_tau,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end subroutine grid_mechanical_spectral_polarisation_updateCoords end subroutine grid_mechanical_spectral_polarization_updateCoords
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file. !> @brief Write current solver and constitutive data for restart to file.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_spectral_polarisation_restartWrite() subroutine grid_mechanical_spectral_polarization_restartWrite()
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
@ -491,7 +497,7 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite()
call DMDAVecRestoreArrayReadF90(DM_mech,FandF_tau_PETSc,FandF_tau,err_PETSc) call DMDAVecRestoreArrayReadF90(DM_mech,FandF_tau_PETSc,FandF_tau,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end subroutine grid_mechanical_spectral_polarisation_restartWrite end subroutine grid_mechanical_spectral_polarization_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -644,4 +650,4 @@ subroutine formResidual(residual_subdomain, FandF_tau, &
end subroutine formResidual end subroutine formResidual
end module grid_mechanical_spectral_polarisation end module grid_mechanical_spectral_polarization

37
src/grid/grid_thermal_spectral.f90
View File

@ -16,6 +16,7 @@ module grid_thermal_spectral
use prec use prec
use parallelization use parallelization
use IO use IO
use misc
use CLI use CLI
use HDF5_utilities use HDF5_utilities
use HDF5 use HDF5
@ -66,9 +67,11 @@ module grid_thermal_spectral
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data !> @brief Allocate all necessary fields and fill them with data, potentially from restart info.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_thermal_spectral_init() subroutine grid_thermal_spectral_init(num_grid)
type(tDict), pointer, intent(in) :: num_grid
integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global
integer :: i, j, k, ce integer :: i, j, k, ce
@ -79,7 +82,11 @@ subroutine grid_thermal_spectral_init()
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
real(pREAL), dimension(1,product(cells(1:2))*cells3) :: tempN real(pREAL), dimension(1,product(cells(1:2))*cells3) :: tempN
type(tDict), pointer :: & type(tDict), pointer :: &
num_grid num_grid_thermal
character(len=:), allocatable :: &
extmsg, &
petsc_options
print'(/,1x,a)', '<<<+- grid_thermal_spectral init -+>>>' print'(/,1x,a)', '<<<+- grid_thermal_spectral init -+>>>'
@ -91,21 +98,23 @@ subroutine grid_thermal_spectral_init()
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict) num_grid_thermal => num_grid%get_dict('thermal',defaultVal=emptyDict)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
num%eps_thermal_atol = num_grid%get_asReal('eps_thermal_atol',defaultVal=1.0e-2_pREAL)
num%eps_thermal_rtol = num_grid%get_asReal('eps_thermal_rtol',defaultVal=1.0e-6_pREAL)
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') num%itmax = num_grid_thermal%get_asInt('N_iter_max', defaultVal=100)
if (num%eps_thermal_atol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_thermal_atol') num%eps_thermal_atol = num_grid_thermal%get_asReal('eps_abs_T', defaultVal=1.0e-2_pREAL)
if (num%eps_thermal_rtol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_thermal_rtol') num%eps_thermal_rtol = num_grid_thermal%get_asReal('eps_rel_T', defaultVal=1.0e-6_pREAL)
extmsg = ''
if (num%eps_thermal_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_T'
if (num%eps_thermal_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_T'
if (num%itmax < 1) extmsg = trim(extmsg)//' N_iter_max'
if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type newtonls -thermal_snes_mf & petsc_options = misc_prefixOptions('-snes_type newtonls -snes_mf -snes_ksp_ew -ksp_type fgmres '// &
&-thermal_snes_ksp_ew -thermal_ksp_type fgmres',err_PETSc) num_grid_thermal%get_asStr('PETSc_options',defaultVal=''), 'thermal_')
CHKERRQ(err_PETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asStr('petsc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

55
src/grid/spectral_utilities.f90
View File

@ -100,12 +100,18 @@ module spectral_utilities
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
DERIVATIVE_CONTINUOUS_ID, & DERIVATIVE_CONTINUOUS_ID, &
DERIVATIVE_CENTRAL_DIFF_ID, & DERIVATIVE_CENTRAL_DIFF_ID, &
DERIVATIVE_FWBW_DIFF_ID DERIVATIVE_FWBW_DIFF_ID, &
DIVERGENCE_CORRECTION_NONE_ID, &
DIVERGENCE_CORRECTION_SIZE_ID, &
DIVERGENCE_CORRECTION_SIZE_GRID_ID
end enum end enum
integer(kind(DERIVATIVE_CONTINUOUS_ID)) :: & integer(kind(DERIVATIVE_CONTINUOUS_ID)) :: &
spectral_derivative_ID spectral_derivative_ID
integer(kind(DIVERGENCE_CORRECTION_NONE_ID)) :: &
divergence_correction_ID
public :: & public :: &
spectral_utilities_init, & spectral_utilities_init, &
utilities_updateGamma, & utilities_updateGamma, &
@ -146,8 +152,9 @@ subroutine spectral_utilities_init()
vectorSize = 3_C_INTPTR_T, & vectorSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T tensorSize = 9_C_INTPTR_T
type(tDict) , pointer :: & type(tDict) , pointer :: &
num_grid num_solver, &
num_grid, &
num_grid_fft
print'(/,1x,a)', '<<<+- spectral_utilities init -+>>>' print'(/,1x,a)', '<<<+- spectral_utilities init -+>>>'
@ -163,8 +170,10 @@ subroutine spectral_utilities_init()
print'( 1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' print'( 1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80' print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict)
num_grid => num_solver%get_dict('grid',defaultVal=emptyDict)
num_grid_fft => num_grid%get_dict('FFT',defaultVal=emptyDict)
num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict)
call PetscOptionsClear(PETSC_NULL_OPTIONS,err_PETSc) call PetscOptionsClear(PETSC_NULL_OPTIONS,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,& call PetscOptionsInsertString(PETSC_NULL_OPTIONS,&
@ -174,13 +183,21 @@ subroutine spectral_utilities_init()
cells1Red = cells(1)/2 + 1 cells1Red = cells(1)/2 + 1
wgt = real(product(cells),pREAL)**(-1) wgt = real(product(cells),pREAL)**(-1)
num%memory_efficient = num_grid%get_asInt('memory_efficient', defaultVal=1) > 0 ! ToDo: should be logical in YAML file num%memory_efficient = num_grid_fft%get_asBool('memory_efficient', defaultVal=.true.)
num%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2)
if (num%divergence_correction < 0 .or. num%divergence_correction > 2) & select case (num_grid_fft%get_asStr('divergence_correction',defaultVal='grid+size'))
call IO_error(301,ext_msg='divergence_correction') case ('none')
divergence_correction_ID = DIVERGENCE_CORRECTION_NONE_ID
case ('size')
divergence_correction_ID = DIVERGENCE_CORRECTION_SIZE_ID
case ('grid+size', 'size+grid')
divergence_correction_ID = DIVERGENCE_CORRECTION_SIZE_GRID_ID
case default
call IO_error(301,ext_msg=trim(num_grid_fft%get_asStr('divergence_correction')))
end select
select case (num_grid%get_asStr('derivative',defaultVal='continuous'))
select case (num_grid_fft%get_asStr('derivative',defaultVal='continuous'))
case ('continuous') case ('continuous')
spectral_derivative_ID = DERIVATIVE_CONTINUOUS_ID spectral_derivative_ID = DERIVATIVE_CONTINUOUS_ID
case ('central_difference') case ('central_difference')
@ -188,18 +205,18 @@ subroutine spectral_utilities_init()
case ('FWBW_difference') case ('FWBW_difference')
spectral_derivative_ID = DERIVATIVE_FWBW_DIFF_ID spectral_derivative_ID = DERIVATIVE_FWBW_DIFF_ID
case default case default
call IO_error(892,ext_msg=trim(num_grid%get_asStr('derivative'))) call IO_error(892,ext_msg=trim(num_grid_fft%get_asStr('derivative')))
end select end select
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and ! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and
! resolution-independent divergence ! resolution-independent divergence
if (num%divergence_correction == 1) then if (divergence_correction_ID == DIVERGENCE_CORRECTION_NONE_ID) then
do j = 1, 3 do j = 1, 3
if (j /= minloc(geomSize,1) .and. j /= maxloc(geomSize,1)) & if (j /= minloc(geomSize,1) .and. j /= maxloc(geomSize,1)) &
scaledGeomSize = geomSize/geomSize(j) scaledGeomSize = geomSize/geomSize(j)
end do end do
elseif (num%divergence_correction == 2) then elseif (divergence_correction_ID == DIVERGENCE_CORRECTION_SIZE_GRID_ID) then
do j = 1, 3 do j = 1, 3
if ( j /= int(minloc(geomSize/real(cells,pREAL),1)) & if ( j /= int(minloc(geomSize/real(cells,pREAL),1)) &
.and. j /= int(maxloc(geomSize/real(cells,pREAL),1))) & .and. j /= int(maxloc(geomSize/real(cells,pREAL),1))) &
@ -209,24 +226,24 @@ subroutine spectral_utilities_init()
scaledGeomSize = geomSize scaledGeomSize = geomSize
end if end if
select case(IO_lc(num_grid%get_asStr('fftw_plan_mode',defaultVal='FFTW_MEASURE'))) select case(IO_lc(num_grid_fft%get_asStr('FFTW_plan_mode',defaultVal='FFTW_MEASURE')))
case('fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution case('fftw_estimate', 'FFTW_ESTIMATE') ! ordered from slow execution (but fast plan creation) to fast execution
FFTW_planner_flag = FFTW_ESTIMATE FFTW_planner_flag = FFTW_ESTIMATE
case('fftw_measure') case('fftw_measure', 'FFTW_MEASURE')
FFTW_planner_flag = FFTW_MEASURE FFTW_planner_flag = FFTW_MEASURE
case('fftw_patient') case('fftw_patient', 'FFTW_PATIENT')
FFTW_planner_flag = FFTW_PATIENT FFTW_planner_flag = FFTW_PATIENT
case('fftw_exhaustive') case('fftw_exhaustive', 'FFTW_EXHAUSTIVE')
FFTW_planner_flag = FFTW_EXHAUSTIVE FFTW_planner_flag = FFTW_EXHAUSTIVE
case default case default
call IO_warning(47,'using default FFTW_MEASURE instead of "'//trim(num_grid%get_asStr('fftw_plan_mode'))//'"') call IO_warning(47,'using default FFTW_MEASURE instead of "'//trim(num_grid_fft%get_asStr('plan_mode'))//'"')
FFTW_planner_flag = FFTW_MEASURE FFTW_planner_flag = FFTW_MEASURE
end select end select
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details) ! general initialization of FFTW (see manual on fftw.org for more details)
if (pREAL /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match' if (pREAL /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match'
call fftw_set_timelimit(num_grid%get_asReal('fftw_timelimit',defaultVal=300.0_pREAL)) call fftw_set_timelimit(num_grid_fft%get_asReal('FFTW_timelimit',defaultVal=300.0_pREAL))
print'(/,1x,a)', 'FFTW initialized'; flush(IO_STDOUT) print'(/,1x,a)', 'FFTW initialized'; flush(IO_STDOUT)

2
src/homogenization.f90
View File

@ -16,7 +16,7 @@ module homogenization
use HDF5 use HDF5
use HDF5_utilities use HDF5_utilities
use result use result
use lattice use crystal
implicit none(type,external) implicit none(type,external)
private private

4
src/homogenization_mechanical_RGC.f90
View File

@ -8,7 +8,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization:mechanical) RGC submodule(homogenization:mechanical) RGC
use rotations use rotations
use lattice use crystal
type :: tParameters type :: tParameters
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -654,7 +654,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
C = phase_homogenizedC66(material_ID_phase(co,ce),material_entry_phase(co,ce)) ! damage not included! C = phase_homogenizedC66(material_ID_phase(co,ce),material_entry_phase(co,ce)) ! damage not included!
equivalentMu = lattice_isotropic_mu(C,'isostrain') equivalentMu = crystal_isotropic_mu(C,'isostrain')
end function equivalentMu end function equivalentMu

4
src/materialpoint.f90
View File

@ -20,7 +20,7 @@ module materialpoint
use rotations use rotations
use polynomials use polynomials
use tables use tables
use lattice use crystal
use material use material
use phase use phase
use homogenization use homogenization
@ -64,7 +64,7 @@ subroutine materialpoint_initAll()
call rotations_init() call rotations_init()
call polynomials_init() call polynomials_init()
call tables_init() call tables_init()
call lattice_init() call crystal_init()
#if defined(MESH) #if defined(MESH)
call discretization_mesh_init(restart=CLI_restartInc>0) call discretization_mesh_init(restart=CLI_restartInc>0)
#elif defined(GRID) #elif defined(GRID)

32
src/misc.f90
View File

@ -5,6 +5,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module misc module misc
use prec use prec
use constants
implicit none(type,external) implicit none(type,external)
private private
@ -18,7 +19,8 @@ module misc
public :: & public :: &
misc_init, & misc_init, &
misc_optional misc_optional, &
misc_prefixOptions
contains contains
@ -110,6 +112,28 @@ pure function misc_optional_str(given,default) result(var)
end function misc_optional_str end function misc_optional_str
!--------------------------------------------------------------------------------------------------
!> @brief Add prefix to options in string.
!> @detail An option starts with a dash followed by at least one letter.
!--------------------------------------------------------------------------------------------------
pure function misc_prefixOptions(string,prefix) result(prefixed)
character(len=*), intent(in) :: string,prefix
character(len=:), allocatable :: prefixed
integer :: i,N
prefixed = ''
N = len(string)
do i = 1, N
prefixed = prefixed//string(i:i)
if (string(i:i) == '-' .and. verify(string(min(i+1,N):min(i+1,N)),LOWER//UPPER) == 0) &
prefixed = prefixed//prefix
end do
end function misc_prefixOptions
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Check correctness of some misc functions. !> @brief Check correctness of some misc functions.
@ -117,6 +141,8 @@ end function misc_optional_str
subroutine misc_selfTest() subroutine misc_selfTest()
real(pREAL) :: r real(pREAL) :: r
character(len=:), allocatable :: str,out
call random_number(r) call random_number(r)
if (test_str('DAMASK') /= 'DAMASK') error stop 'optional_str, present' if (test_str('DAMASK') /= 'DAMASK') error stop 'optional_str, present'
@ -132,6 +158,10 @@ subroutine misc_selfTest()
if (.not. test_bool()) error stop 'optional_bool, not present' if (.not. test_bool()) error stop 'optional_bool, not present'
if (misc_optional(default=r>0.5_pREAL) .neqv. r>0.5_pREAL) error stop 'optional_bool, default only' if (misc_optional(default=r>0.5_pREAL) .neqv. r>0.5_pREAL) error stop 'optional_bool, default only'
str='-a -1 -more 123 -flag -'
out=misc_prefixOptions(str,'p_')
if (out /= '-p_a -1 -p_more 123 -p_flag -') error stop 'misc_prefixOptions'
contains contains
function test_str(str_in) result(str_out) function test_str(str_in) result(str_out)

4
src/phase.f90
View File

@ -14,7 +14,7 @@ module phase
use config use config
use material use material
use result use result
use lattice use crystal
use discretization use discretization
use parallelization use parallelization
use HDF5 use HDF5
@ -358,7 +358,7 @@ module phase
config, & config, &
material, & material, &
result, & result, &
lattice, & crystal, &
discretization, & discretization, &
HDF5_utilities HDF5_utilities
#endif #endif

2
src/phase_damage_anisobrittle.f90
View File

@ -77,7 +77,7 @@ module function anisobrittle_init() result(mySources)
prm%s_crit = src%get_as1dReal('s_crit',requiredSize=size(N_cl)) prm%s_crit = src%get_as1dReal('s_crit',requiredSize=size(N_cl))
prm%g_crit = src%get_as1dReal('g_crit',requiredSize=size(N_cl)) prm%g_crit = src%get_as1dReal('g_crit',requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase_lattice(ph),phase_cOverA(ph)) prm%cleavage_systems = crystal_SchmidMatrix_cleavage(N_cl,phase_lattice(ph),phase_cOverA(ph))
! expand: family => system ! expand: family => system
prm%s_crit = math_expand(prm%s_crit,N_cl) prm%s_crit = math_expand(prm%s_crit,N_cl)

2
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -92,7 +92,7 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
Alpha = 0.0_pREAL Alpha = 0.0_pREAL
Alpha(1,1) = prm%Alpha_11%at(T) Alpha(1,1) = prm%Alpha_11%at(T)
if (any(phase_lattice(ph) == ['hP','tI'])) Alpha(3,3) = prm%Alpha_33%at(T) if (any(phase_lattice(ph) == ['hP','tI'])) Alpha(3,3) = prm%Alpha_33%at(T)
Alpha = lattice_symmetrize_33(Alpha,phase_lattice(ph)) Alpha = crystal_symmetrize_33(Alpha,phase_lattice(ph))
Li = dot_T * Alpha Li = dot_T * Alpha
end associate end associate

6
src/phase_mechanical_elastic.f90
View File

@ -97,7 +97,7 @@ pure module function elastic_C66(ph,en) result(C66)
if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66%at(T) if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66%at(T)
C66 = lattice_symmetrize_C66(C66,phase_lattice(ph)) C66 = crystal_symmetrize_C66(C66,phase_lattice(ph))
end associate end associate
@ -119,7 +119,7 @@ pure module function elastic_mu(ph,en,isotropic_bound) result(mu)
associate(prm => param(ph)) associate(prm => param(ph))
mu = lattice_isotropic_mu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph)) mu = crystal_isotropic_mu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph))
end associate end associate
@ -141,7 +141,7 @@ pure module function elastic_nu(ph,en,isotropic_bound) result(nu)
associate(prm => param(ph)) associate(prm => param(ph))
nu = lattice_isotropic_nu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph)) nu = crystal_isotropic_nu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph))
end associate end associate

14
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -149,13 +149,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%P_sl = crystal_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
if (phase_lattice(ph) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dReal('a_nonSchmid',defaultVal = emptyRealArray) a = pl%get_as1dReal('a_nonSchmid',defaultVal = emptyRealArray)
prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = crystal_nonSchmidMatrix(N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = crystal_nonSchmidMatrix(N_sl,a,-1)
else else
prm%P_nS_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
@ -184,13 +184,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%d_caron = prm%b_sl * pl%get_asReal('D_a') prm%d_caron = prm%b_sl * pl%get_asReal('D_a')
prm%f_at = prm%b_sl**3*pl%get_asReal('f_at') prm%f_at = prm%b_sl**3*pl%get_asReal('f_at')
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), & prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), &
phase_lattice(ph)) phase_lattice(ph))
prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) & prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) & * crystal_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) &
+ spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) & + spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph)) * crystal_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph))
! sanity checks ! sanity checks
if ( prm%D_0 < 0.0_pREAL) extmsg = trim(extmsg)//' D_0' if ( prm%D_0 < 0.0_pREAL) extmsg = trim(extmsg)//' D_0'

36
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -73,7 +73,7 @@ submodule(phase:plastic) dislotwin
integer, allocatable, dimension(:,:) :: & integer, allocatable, dimension(:,:) :: &
fcc_twinNucleationSlipPair ! ToDo: Better name? Is also used for trans fcc_twinNucleationSlipPair ! ToDo: Better name? Is also used for trans
character(len=:), allocatable :: & character(len=:), allocatable :: &
lattice_tr, & crystal_tr, &
isotropic_bound isotropic_bound
character(len=pSTRLEN), allocatable, dimension(:) :: & character(len=pSTRLEN), allocatable, dimension(:) :: &
output output
@ -202,9 +202,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%P_sl = crystal_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%n0_sl = lattice_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%n0_sl = crystal_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%extendedDislocations = pl%get_asBool('extend_dislocations',defaultVal=.false.) prm%extendedDislocations = pl%get_asBool('extend_dislocations',defaultVal=.false.)
prm%omitDipoles = pl%get_asBool('omit_dipoles', defaultVal=.false.) prm%omitDipoles = pl%get_asBool('omit_dipoles', defaultVal=.false.)
@ -226,15 +226,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl) defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl)
prm%d_caron = prm%b_sl * pl%get_asReal('D_a') prm%d_caron = prm%b_sl * pl%get_asReal('D_a')
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph)) prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph))
prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) & prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) & * crystal_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) &
+ spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) & + spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph)) * crystal_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. N_sl(1) == 12 prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. N_sl(1) == 12
if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_CF_TWINNUCLEATIONSLIPPAIR if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = crystal_CF_TWINNUCLEATIONSLIPPAIR
! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex)
! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 ! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
@ -274,9 +274,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray) prm%N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(prm%N_tw)) prm%sum_N_tw = sum(abs(prm%N_tw))
twinActive: if (prm%sum_N_tw > 0) then twinActive: if (prm%sum_N_tw > 0) then
prm%systems_tw = lattice_labels_twin(prm%N_tw,phase_lattice(ph)) prm%systems_tw = crystal_labels_twin(prm%N_tw,phase_lattice(ph))
prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%P_tw = crystal_SchmidMatrix_twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%gamma_char_tw = lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%gamma_char_tw = crystal_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%L_tw = pl%get_asReal('L_tw') prm%L_tw = pl%get_asReal('L_tw')
prm%i_tw = pl%get_asReal('i_tw') prm%i_tw = pl%get_asReal('i_tw')
@ -285,7 +285,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%t_tw = math_expand(pl%get_as1dReal('t_tw', requiredSize=size(prm%N_tw)),prm%N_tw) prm%t_tw = math_expand(pl%get_as1dReal('t_tw', requiredSize=size(prm%N_tw)),prm%N_tw)
prm%r = math_expand(pl%get_as1dReal('p_tw', requiredSize=size(prm%N_tw)),prm%N_tw) prm%r = math_expand(pl%get_as1dReal('p_tw', requiredSize=size(prm%N_tw)),prm%N_tw)
prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,pl%get_as1dReal('h_tw-tw'), & prm%h_tw_tw = crystal_interaction_TwinByTwin(prm%N_tw,pl%get_as1dReal('h_tw-tw'), &
phase_lattice(ph)) phase_lattice(ph))
! sanity checks ! sanity checks
@ -309,7 +309,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray) prm%N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray)
prm%sum_N_tr = sum(abs(prm%N_tr)) prm%sum_N_tr = sum(abs(prm%N_tr))
transActive: if (prm%sum_N_tr > 0) then transActive: if (prm%sum_N_tr > 0) then
prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr,'hP',prm%cOverA_hP) prm%P_tr = crystal_SchmidMatrix_trans(prm%N_tr,'hP',prm%cOverA_hP)
prm%Delta_G = polynomial(pl,'Delta_G','T') prm%Delta_G = polynomial(pl,'Delta_G','T')
prm%i_tr = pl%get_asReal('i_tr') prm%i_tr = pl%get_asReal('i_tr')
@ -324,7 +324,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
a_cF = prm%b_tr(1)*sqrt(6.0_pREAL) ! b_tr is Shockley partial a_cF = prm%b_tr(1)*sqrt(6.0_pREAL) ! b_tr is Shockley partial
prm%h = 5.0_pREAL * a_cF/sqrt(3.0_pREAL) prm%h = 5.0_pREAL * a_cF/sqrt(3.0_pREAL)
prm%rho = 4.0_pREAL/(sqrt(3.0_pREAL)*a_cF**2)/N_A prm%rho = 4.0_pREAL/(sqrt(3.0_pREAL)*a_cF**2)/N_A
prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dReal('h_tr-tr'),& prm%h_tr_tr = crystal_interaction_TransByTrans(prm%N_tr,pl%get_as1dReal('h_tr-tr'),&
phase_lattice(ph)) phase_lattice(ph))
@ -372,13 +372,13 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%Gamma_sf = polynomial(pl,'Gamma_sf','T') prm%Gamma_sf = polynomial(pl,'Gamma_sf','T')
slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dReal('h_sl-tw'), & prm%h_sl_tw = crystal_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dReal('h_sl-tw'), &
phase_lattice(ph)) phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation' if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
end if slipAndTwinActive end if slipAndTwinActive
slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dReal('h_sl-tr'), & prm%h_sl_tr = crystal_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dReal('h_sl-tr'), &
phase_lattice(ph)) phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation' if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
end if slipAndTransActive end if slipAndTransActive
@ -480,7 +480,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
homogenizedC = f_matrix * C homogenizedC = f_matrix * C
twinActive: if (prm%sum_N_tw > 0) then twinActive: if (prm%sum_N_tw > 0) then
C66_tw = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph)) C66_tw = crystal_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph))
do i = 1, prm%sum_N_tw do i = 1, prm%sum_N_tw
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tw(i,en)*C66_tw(1:6,1:6,i) + stt%f_tw(i,en)*C66_tw(1:6,1:6,i)
@ -488,7 +488,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
end if twinActive end if twinActive
transActive: if (prm%sum_N_tr > 0) then transActive: if (prm%sum_N_tr > 0) then
C66_tr = lattice_C66_trans(prm%N_tr,C,'hP',prm%cOverA_hP) C66_tr = crystal_C66_trans(prm%N_tr,C,'hP',prm%cOverA_hP)
do i = 1, prm%sum_N_tr do i = 1, prm%sum_N_tr
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tr(i,en)*C66_tr(1:6,1:6,i) + stt%f_tr(i,en)*C66_tr(1:6,1:6,i)

10
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -139,14 +139,14 @@ module function plastic_kinehardening_init() result(myPlasticity)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
prm%P = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%P = crystal_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
if (phase_lattice(ph) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dReal('a_nonSchmid',defaultVal=emptyRealArray) a = pl%get_as1dReal('a_nonSchmid',defaultVal=emptyRealArray)
prm%nonSchmidActive = size(a) > 0 prm%nonSchmidActive = size(a) > 0
prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = crystal_nonSchmidMatrix(N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = crystal_nonSchmidMatrix(N_sl,a,-1)
else else
prm%P_nS_pos = prm%P prm%P_nS_pos = prm%P
prm%P_nS_neg = prm%P prm%P_nS_neg = prm%P
@ -155,7 +155,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
prm%dot_gamma_0 = pl%get_asReal('dot_gamma_0') prm%dot_gamma_0 = pl%get_asReal('dot_gamma_0')
prm%n = pl%get_asReal('n') prm%n = pl%get_asReal('n')
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), & prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), &
phase_lattice(ph)) phase_lattice(ph))
xi_0 = math_expand(pl%get_as1dReal('xi_0', requiredSize=size(N_sl)),N_sl) xi_0 = math_expand(pl%get_as1dReal('xi_0', requiredSize=size(N_sl)),N_sl)

20
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -249,30 +249,30 @@ module function plastic_nonlocal_init() result(myPlasticity)
ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(ini%N_sl)) prm%sum_N_sl = sum(abs(ini%N_sl))
slipActive: if (prm%sum_N_sl > 0) then slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(ini%N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(ini%N_sl,phase_lattice(ph))
prm%P_sl = lattice_SchmidMatrix_slip(ini%N_sl,phase_lattice(ph), phase_cOverA(ph)) prm%P_sl = crystal_SchmidMatrix_slip(ini%N_sl,phase_lattice(ph), phase_cOverA(ph))
if (phase_lattice(ph) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dReal('a_nonSchmid',defaultVal = emptyRealArray) a = pl%get_as1dReal('a_nonSchmid',defaultVal = emptyRealArray)
if (size(a) > 0) prm%nonSchmidActive = .true. if (size(a) > 0) prm%nonSchmidActive = .true.
prm%P_nS_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1) prm%P_nS_pos = crystal_nonSchmidMatrix(ini%N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1) prm%P_nS_neg = crystal_nonSchmidMatrix(ini%N_sl,a,-1)
else else
prm%P_nS_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
end if end if
prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dReal('h_sl-sl'), & prm%h_sl_sl = crystal_interaction_SlipBySlip(ini%N_sl,pl%get_as1dReal('h_sl-sl'), &
phase_lattice(ph)) phase_lattice(ph))
prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase_lattice(ph),& prm%forestProjection_edge = crystal_forestProjection_edge (ini%N_sl,phase_lattice(ph),&
phase_cOverA(ph)) phase_cOverA(ph))
prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase_lattice(ph),& prm%forestProjection_screw = crystal_forestProjection_screw(ini%N_sl,phase_lattice(ph),&
phase_cOverA(ph)) phase_cOverA(ph))
prm%slip_direction = lattice_slip_direction (ini%N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%slip_direction = crystal_slip_direction (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%slip_transverse = crystal_slip_transverse(ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%slip_normal = lattice_slip_normal (ini%N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%slip_normal = crystal_slip_normal (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
! collinear systems (only for octahedral slip systems in fcc) ! collinear systems (only for octahedral slip systems in fcc)
allocate(prm%colinearSystem(prm%sum_N_sl), source = -1) allocate(prm%colinearSystem(prm%sum_N_sl), source = -1)

22
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -149,21 +149,21 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
prm%h_int = math_expand(pl%get_as1dReal('h_int', requiredSize=size(N_sl), & prm%h_int = math_expand(pl%get_as1dReal('h_int', requiredSize=size(N_sl), &
defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl) defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl)
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph)) prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) prm%P_sl = crystal_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
if (phase_lattice(ph) == 'cI') then if (phase_lattice(ph) == 'cI') then
a = pl%get_as1dReal('a_nonSchmid',defaultVal=emptyRealArray) a = pl%get_as1dReal('a_nonSchmid',defaultVal=emptyRealArray)
if (size(a) > 0) prm%nonSchmidActive = .true. if (size(a) > 0) prm%nonSchmidActive = .true.
prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%P_nS_pos = crystal_nonSchmidMatrix(N_sl,a,+1)
prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%P_nS_neg = crystal_nonSchmidMatrix(N_sl,a,-1)
else else
prm%P_nS_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
end if end if
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
! sanity checks ! sanity checks
if ( prm%dot_gamma_0_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_sl' if ( prm%dot_gamma_0_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_sl'
@ -196,11 +196,11 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
xi_0_tw = math_expand(pl%get_as1dReal('xi_0_tw',requiredSize=size(N_tw)),N_tw) xi_0_tw = math_expand(pl%get_as1dReal('xi_0_tw',requiredSize=size(N_tw)),N_tw)
prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%gamma_char = crystal_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph)) prm%h_tw_tw = crystal_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph))
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%P_tw = crystal_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%systems_tw = lattice_labels_twin(N_tw,phase_lattice(ph)) prm%systems_tw = crystal_labels_twin(N_tw,phase_lattice(ph))
! sanity checks ! sanity checks
if (prm%dot_gamma_0_tw <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_tw' if (prm%dot_gamma_0_tw <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_tw'
@ -216,9 +216,9 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
! slip-twin related parameters ! slip-twin related parameters
slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h_0_tw_sl = pl%get_asReal('h_0_tw-sl') prm%h_0_tw_sl = pl%get_asReal('h_0_tw-sl')
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dReal('h_sl-tw'), & prm%h_sl_tw = crystal_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dReal('h_sl-tw'), &
phase_lattice(ph)) phase_lattice(ph))
prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dReal('h_tw-sl'), & prm%h_tw_sl = crystal_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dReal('h_tw-sl'), &
phase_lattice(ph)) phase_lattice(ph))
else slipAndTwinActive else slipAndTwinActive
allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0 allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0

2
src/phase_thermal.f90
View File

@ -104,7 +104,7 @@ module subroutine thermal_init(phases)
param(ph)%C_p = thermal%get_asReal('C_p') param(ph)%C_p = thermal%get_asReal('C_p')
param(ph)%K(1,1) = thermal%get_asReal('K_11') param(ph)%K(1,1) = thermal%get_asReal('K_11')
if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asReal('K_33') if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asReal('K_33')
param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph)) param(ph)%K = crystal_symmetrize_33(param(ph)%K,phase_lattice(ph))
#if defined(__GFORTRAN__) #if defined(__GFORTRAN__)
param(ph)%output = output_as1dStr(thermal) param(ph)%output = output_as1dStr(thermal)