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

added new, flexible debugging scheme. 

now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing
This commit is contained in:
Martin Diehl 2012-03-08 20:25:28 +00:00
parent dec9451b1e
commit bd9667bd4b
20 changed files with 2790 additions and 2727 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

215
code/CPFEM.f90
View File

@ -23,8 +23,7 @@ MODULE CPFEM
!##############################################################
! *** CPFEM engine ***
!
use prec, only: pReal, &
pInt
use prec, only: pReal
implicit none
real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, &
@ -47,8 +46,8 @@ CONTAINS
subroutine CPFEM_initAll(Temperature,element,IP)
use prec, only: pReal, &
prec_init
use prec, only: prec_init, &
pInt
use numerics, only: numerics_init
use debug, only: debug_init
use FEsolving, only: FE_init
@ -61,8 +60,8 @@ subroutine CPFEM_initAll(Temperature,element,IP)
use homogenization, only: homogenization_init
use IO, only: IO_init
use DAMASK_interface
implicit none
implicit none
integer(pInt), intent(in) :: element, & ! FE element number
IP ! FE integration point number
real(pReal), intent(in) :: Temperature ! temperature
@ -79,19 +78,19 @@ subroutine CPFEM_initAll(Temperature,element,IP)
n = n+1_pInt
if (.not. CPFEM_init_inProgress) then ! yes my thread won!
CPFEM_init_inProgress = .true.
call prec_init()
call IO_init()
call numerics_init()
call debug_init()
call math_init()
call FE_init()
call prec_init
call IO_init
call numerics_init
call debug_init
call math_init
call FE_init
call mesh_init(IP, element) ! pass on coordinates to alter calcMode of first ip
call lattice_init()
call material_init()
call constitutive_init()
call lattice_init
call material_init
call constitutive_init
call crystallite_init(Temperature) ! (have to) use temperature of first IP for whole model
call homogenization_init(Temperature)
call CPFEM_init()
call CPFEM_init
if (trim(FEsolver)/='Spectral') call DAMASK_interface_init() ! Spectral solver is doing initialization earlier
CPFEM_init_done = .true.
CPFEM_init_inProgress = .false.
@ -101,18 +100,20 @@ subroutine CPFEM_initAll(Temperature,element,IP)
endif
endif
end subroutine
end subroutine CPFEM_initAll
!*********************************************************
!*** allocate the arrays defined in module CPFEM ***
!*** and initialize them ***
!*********************************************************
subroutine CPFEM_init()
subroutine CPFEM_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt
use debug, only: debug_verbosity
use debug, only: debug_what, &
debug_CPFEM, &
debug_levelBasic
use IO, only: IO_read_jobBinaryFile
use FEsolving, only: parallelExecution, &
symmetricSolver, &
@ -133,7 +134,6 @@ subroutine CPFEM_init()
implicit none
integer(pInt) i,j,k,l,m
! initialize stress and jacobian to zero
@ -143,7 +143,7 @@ subroutine CPFEM_init()
! *** restore the last converged values of each essential variable from the binary file
if (restartRead) then
if (debug_verbosity > 0) then
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a)') '<< CPFEM >> Restored state variables of last converged step from binary files'
!$OMP END CRITICAL (write2out)
@ -205,7 +205,7 @@ subroutine CPFEM_init()
write(6,*) '<<<+- cpfem init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
if (debug_verbosity > 0) then
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0) then
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs)
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
@ -216,7 +216,7 @@ subroutine CPFEM_init()
call flush(6)
!$OMP END CRITICAL (write2out)
endsubroutine
end subroutine CPFEM_init
!***********************************************************************
@ -228,81 +228,82 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE
!*** variables and functions from other modules ***!
use prec, only: pReal, &
pInt
use numerics, only: defgradTolerance, &
iJacoStiffness
use debug, only: debug_e, &
debug_i, &
debug_selectiveDebugger, &
debug_verbosity, &
debug_stressMaxLocation, &
debug_stressMinLocation, &
debug_jacobianMaxLocation, &
debug_jacobianMinLocation, &
debug_stressMax, &
debug_stressMin, &
debug_jacobianMax, &
debug_jacobianMin
use FEsolving, only: parallelExecution, &
outdatedFFN1, &
terminallyIll, &
cycleCounter, &
theInc, &
theTime, &
theDelta, &
FEsolving_execElem, &
FEsolving_execIP, &
restartWrite
use math, only: math_identity2nd, &
math_mul33x33, &
math_det33, &
math_transpose33, &
math_I3, &
math_Mandel3333to66, &
math_Mandel66to3333, &
math_Mandel33to6, &
math_Mandel6to33
use mesh, only: mesh_FEasCP, &
mesh_NcpElems, &
mesh_maxNips, &
mesh_element, &
mesh_node0, &
mesh_node, &
mesh_ipCenterOfGravity, &
mesh_build_subNodeCoords, &
mesh_build_ipVolumes, &
mesh_build_ipCoordinates, &
FE_Nips, &
FE_Nnodes
use material, only: homogenization_maxNgrains, &
microstructure_elemhomo, &
material_phase
use constitutive, only: constitutive_state0,constitutive_state
use crystallite, only: crystallite_partionedF,&
crystallite_F0, &
crystallite_Fp0, &
crystallite_Fp, &
crystallite_Lp0, &
crystallite_Lp, &
crystallite_dPdF0, &
crystallite_dPdF, &
crystallite_Tstar0_v, &
crystallite_Tstar_v
use homogenization, only: homogenization_sizeState, &
homogenization_state, &
homogenization_state0, &
materialpoint_F, &
materialpoint_F0, &
materialpoint_P, &
materialpoint_dPdF, &
materialpoint_results, &
materialpoint_sizeResults, &
materialpoint_Temperature, &
materialpoint_stressAndItsTangent, &
materialpoint_postResults
use IO, only: IO_write_jobBinaryFile, &
IO_warning
use prec, only: pInt
use numerics, only: defgradTolerance, &
iJacoStiffness
use debug, only: debug_what, &
debug_CPFEM, &
debug_levelBasic, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_stressMaxLocation, &
debug_stressMinLocation, &
debug_jacobianMaxLocation, &
debug_jacobianMinLocation, &
debug_stressMax, &
debug_stressMin, &
debug_jacobianMax, &
debug_jacobianMin
use FEsolving, only: parallelExecution, &
outdatedFFN1, &
terminallyIll, &
cycleCounter, &
theInc, &
theTime, &
theDelta, &
FEsolving_execElem, &
FEsolving_execIP, &
restartWrite
use math, only: math_identity2nd, &
math_mul33x33, &
math_det33, &
math_transpose33, &
math_I3, &
math_Mandel3333to66, &
math_Mandel66to3333, &
math_Mandel33to6, &
math_Mandel6to33
use mesh, only: mesh_FEasCP, &
mesh_NcpElems, &
mesh_maxNips, &
mesh_element, &
mesh_node0, &
mesh_node, &
mesh_ipCenterOfGravity, &
mesh_build_subNodeCoords, &
mesh_build_ipVolumes, &
mesh_build_ipCoordinates, &
FE_Nips, &
FE_Nnodes
use material, only: homogenization_maxNgrains, &
microstructure_elemhomo, &
material_phase
use constitutive, only: constitutive_state0,constitutive_state
use crystallite, only: crystallite_partionedF,&
crystallite_F0, &
crystallite_Fp0, &
crystallite_Fp, &
crystallite_Lp0, &
crystallite_Lp, &
crystallite_dPdF0, &
crystallite_dPdF, &
crystallite_Tstar0_v, &
crystallite_Tstar_v
use homogenization, only: homogenization_sizeState, &
homogenization_state, &
homogenization_state0, &
materialpoint_F, &
materialpoint_F0, &
materialpoint_P, &
materialpoint_dPdF, &
materialpoint_results, &
materialpoint_sizeResults, &
materialpoint_Temperature, &
materialpoint_stressAndItsTangent, &
materialpoint_postResults
use IO, only: IO_write_jobBinaryFile, &
IO_warning
use DAMASK_interface
implicit none
@ -359,7 +360,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
cp_en = mesh_FEasCP('elem',element)
if (debug_verbosity > 0 .and. cp_en == 1 .and. IP == 1) then
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt .and. cp_en == 1 .and. IP == 1) then
!$OMP CRITICAL (write2out)
write(6,*)
write(6,'(a)') '#############################################'
@ -396,7 +397,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
j = 1:mesh_maxNips, &
k = 1:mesh_NcpElems ) &
constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites
if (debug_verbosity > 0) then
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a)') '<< CPFEM >> Aging states'
if (debug_e == cp_en .and. debug_i == IP) then
@ -418,7 +419,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
! * dump the last converged values of each essential variable to a binary file
if (restartWrite) then
if (debug_verbosity > 0) then
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a)') '<< CPFEM >> Writing state variables of last converged step to binary files'
!$OMP END CRITICAL (write2out)
@ -487,7 +488,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then
if (.not. terminallyIll .and. .not. outdatedFFN1) then
if (debug_verbosity > 0) then
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at element ip',cp_en,IP
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en))
@ -514,7 +515,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
FEsolving_execElem(2) = cp_en
FEsolving_execIP(1,cp_en) = IP
FEsolving_execIP(2,cp_en) = IP
if (debug_verbosity > 0) then
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,1x,i2)') '<< CPFEM >> Calculation for element ip ',cp_en,IP
!$OMP END CRITICAL (write2out)
@ -525,7 +526,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
!* parallel computation and calulation not yet done
elseif (.not. CPFEM_calc_done) then
if (debug_verbosity > 0) then
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
!$OMP END CRITICAL (write2out)
@ -534,7 +535,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
call mesh_build_subNodeCoords() ! update subnodal coordinates
call mesh_build_ipCoordinates() ! update ip coordinates
endif
if (debug_verbosity > 0) then
if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,a,i8)') '<< CPFEM >> Start stress and tangent ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
!$OMP END CRITICAL (write2out)
@ -640,7 +641,9 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition.
endif
if (mode < 3 .and. debug_verbosity > 0 .and. ((debug_e == cp_en .and. debug_i == IP) .or. .not. debug_selectiveDebugger)) then
if (mode < 3 .and. iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt &
.and. ((debug_e == cp_en .and. debug_i == IP) &
.or. .not. iand(debug_what(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') '<< CPFEM >> stress/MPa at el ip ', cp_en, IP, cauchyStress/1.0e6_pReal
write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') '<< CPFEM >> jacobian/GPa at el ip ', cp_en, IP, transpose(jacobian)/1.0e9_pReal
@ -679,6 +682,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
endif
endif
end subroutine
end subroutine CPFEM_general
END MODULE CPFEM

51
code/DAMASK_marc.f90
View File

@ -58,16 +58,17 @@
#include "prec.f90"
MODULE DAMASK_interface
module DAMASK_interface
character(len=64), parameter :: FEsolver = 'Marc'
character(len=4), parameter :: InputFileExtension = '.dat'
character(len=4), parameter :: LogFileExtension = '.log'
character(len=64), parameter :: FEsolver = 'Marc'
character(len=4), parameter :: InputFileExtension = '.dat'
character(len=4), parameter :: LogFileExtension = '.log'
CONTAINS
contains
subroutine DAMASK_interface_init()
subroutine DAMASK_interface_init
implicit none
!$OMP CRITICAL (write2out)
write(6,*)
@ -75,11 +76,14 @@ subroutine DAMASK_interface_init()
write(6,*) '$Id$'
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
return
end subroutine
end subroutine DAMASK_interface_init
function getSolverWorkingDirectoryName()
implicit none
character(1024) getSolverWorkingDirectoryName, outName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
@ -88,26 +92,27 @@ function getSolverWorkingDirectoryName()
inquire(6, name=outName) ! determine outputfile
getSolverWorkingDirectoryName=outName(1:scan(outName,pathSep,back=.true.))
! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName
end function
end function getSolverWorkingDirectoryName
function getModelName()
use prec
implicit none
character(1024) getModelName
character(1024) :: getModelName
getModelName = getSolverJobName()
end function
end function getModelName
function getSolverJobName()
use prec
use prec, only: pInt
implicit none
character(1024) getSolverJobName, outName
character(1024) :: getSolverJobName, outName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
integer(pInt) extPos
integer(pInt) :: extPos
getSolverJobName=''
outName=''
@ -115,9 +120,10 @@ function getSolverJobName()
extPos = len_trim(outName)-4
getSolverJobName=outName(scan(outName,pathSep,back=.true.)+1:extPos)
! write(6,*) 'getSolverJobName', getSolverJobName
end function
END MODULE
end function getSolverJobName
end module DAMASK_interface
#include "IO.f90"
#include "numerics.f90"
@ -234,8 +240,8 @@ subroutine hypela2(&
use CPFEM, only: CPFEM_initAll,CPFEM_general,CPFEM_init_done
!$ use OMP_LIB ! the openMP function library
!$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS
implicit none
implicit none
! ** Start of generated type statements **
real(pReal) coord, d, de, disp, dispt, dt, e, eigvn, eigvn1, ffn, ffn1
real(pReal) frotn, frotn1, g
@ -365,7 +371,7 @@ subroutine hypela2(&
!$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value
end subroutine
end subroutine hypela2
!********************************************************************
@ -394,8 +400,8 @@ subroutine plotv(&
use mesh, only: mesh_FEasCP
use IO, only: IO_error
use homogenization, only: materialpoint_results,materialpoint_sizeResults
implicit none
implicit none
real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*)
real(pReal) v, t(*)
integer(pInt) m, nn, layer, ndi, nshear, jpltcd
@ -403,6 +409,5 @@ subroutine plotv(&
if (jpltcd > materialpoint_sizeResults) call IO_error(700_pInt,jpltcd) ! complain about out of bounds error
v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m))
return
end subroutine
end subroutine plotv

30
code/DAMASK_spectral_interface.f90
View File

@ -27,13 +27,13 @@ module DAMASK_interface
implicit none
private
character(len=64), parameter, public :: FEsolver = 'Spectral' !> Keyword for spectral solver
character(len=5), parameter, public :: inputFileExtension = '.geom' !> File extension for geometry description
character(len=4), parameter, public :: logFileExtension = '.log' !> Dummy variable as the spectral solver has no log
character(len=1024), private :: geometryParameter, &
loadcaseParameter
character(len=64), parameter, public :: FEsolver = 'Spectral' !< Keyword for spectral solver
character(len=5), parameter, public :: inputFileExtension = '.geom' !< File extension for geometry description
character(len=4), parameter, public :: logFileExtension = '.log' !< Dummy variable as the spectral solver has no log
character(len=1024), private :: geometryParameter, & !< Interpretated parameter given at command line
loadcaseParameter !< Interpretated parameter given at command line
public :: getSolverWorkingDirectoryName, &
public :: getSolverWorkingDirectoryName, & !< Interpretated parameter given at command line
getSolverJobName, &
getLoadCase, &
getLoadCaseName, &
@ -46,7 +46,7 @@ module DAMASK_interface
contains
!--------------------------------------------------------------------------------------------------
!> @brief Initializes the solver by interpreting the command line arguments. Also writes
!> @brief initializes the solver by interpreting the command line arguments. Also writes
!! information on computation on screen
!--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init
@ -54,12 +54,12 @@ subroutine DAMASK_interface_init
use prec, only: pInt
implicit none
character(len=1024) :: commandLine, & !> command line call as string
hostName, & !> name of computer
userName !> name of user calling the executable
character(len=1024) :: commandLine, & !< command line call as string
hostName, & !< name of computer
userName !< name of user calling the executable
integer :: i, &
start = 0,&
length=0
start ,&
length
integer, dimension(8) :: dateAndTime ! type default integer
call get_command(commandLine)
@ -186,11 +186,11 @@ end subroutine DAMASK_interface_init
!--------------------------------------------------------------------------------------------------
!> @brief extract working directory from loadcase file possibly based on current working dir
!--------------------------------------------------------------------------------------------------
function getSolverWorkingDirectoryName()
character(len=1024) function getSolverWorkingDirectoryName()
implicit none
character(len=1024) :: cwd, getSolverWorkingDirectoryName
character :: pathSep
character(len=1024) :: cwd
character :: pathSep
pathSep = getPathSep()

85
code/FEsolving.f90
View File

@ -19,50 +19,81 @@
!##############################################################
!* $Id$
!##############################################################
MODULE FEsolving
module FEsolving
!##############################################################
use prec, only: pInt,pReal
implicit none
integer(pInt) :: &
cycleCounter = 0_pInt, &
theInc = -1_pInt, &
restartInc = 1_pInt
integer(pInt) :: cycleCounter = 0_pInt, theInc = -1_pInt, restartInc = 1_pInt
real(pReal) :: theTime = 0.0_pReal, theDelta = 0.0_pReal
logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.,terminallyIll = .false.
logical :: symmetricSolver = .false.
logical :: parallelExecution = .true.
logical :: restartWrite = .false.
logical :: restartRead = .false.
logical :: lastMode = .true., cutBack = .false.
logical, dimension(:,:), allocatable :: calcMode
integer(pInt), dimension(:,:), allocatable :: FEsolving_execIP
integer(pInt), dimension(2) :: FEsolving_execElem
character(len=1024) FEmodelGeometry
real(pReal) :: &
theTime = 0.0_pReal, &
theDelta = 0.0_pReal
CONTAINS
logical :: &
lastIncConverged = .false., &
outdatedByNewInc = .false., &
outdatedFFN1 = .false., &
terminallyIll = .false., &
symmetricSolver = .false., &
parallelExecution = .true., &
restartWrite = .false., &
restartRead = .false., &
lastMode = .true., &
cutBack = .false.
integer(pInt), dimension(:,:), allocatable :: &
FEsolving_execIP
integer(pInt), dimension(2) :: &
FEsolving_execElem
character(len=1024) :: &
FEmodelGeometry
logical, dimension(:,:), allocatable :: &
calcMode
public :: FE_init
contains
!***********************************************************
! determine whether a symmetric solver is used
! and whether restart is requested
!***********************************************************
subroutine FE_init()
subroutine FE_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt
use debug, only: debug_verbosity
use debug, only: debug_what, &
debug_FEsolving, &
debug_levelBasic
use IO, only: IO_open_inputFile, &
IO_stringPos, &
IO_stringValue, &
IO_intValue, &
IO_lc, &
IO_open_logFile, &
IO_warning
use DAMASK_interface
use IO
implicit none
implicit none
integer(pInt), parameter :: fileunit = 222_pInt
integer(pInt), parameter :: maxNchunks = 6_pInt
integer :: i, start = 0, length=0
integer :: i, start = 0, length ! is save for FE_init (only called once)
integer(pInt) :: j
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
character(len=64) tag
character(len=1024) line, commandLine
character(len=64) :: tag
character(len=1024) :: line, &
commandLine
FEmodelGeometry = getModelName()
call IO_open_inputFile(fileunit,FEmodelGeometry)
if (trim(FEsolver) == 'Spectral') then
call get_command(commandLine) ! may contain uppercase
do i=1,len(commandLine)
@ -73,7 +104,7 @@
start = index(commandLine,'-r ',.true.) + 3 ! set to position after trailing space
if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart
start = index(commandLine,'--restart ',.true.) + 10 ! set to position after trailing space
if(start /= 0_pInt) then ! found something
if(start /= 0) then ! found something
length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument?
read(commandLine(start:start+length),'(I12)') restartInc ! read argument
restartRead = restartInc > 0_pInt
@ -148,7 +179,7 @@
write(6,*) '<<<+- FEsolving init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
if (debug_verbosity > 0) then
if (iand(debug_what(debug_FEsolving),debug_levelBasic) /= 0_pInt) then
write(6,*) 'restart writing: ', restartWrite
write(6,*) 'restart reading: ', restartRead
if (restartRead) write(6,*) 'restart Job: ', trim(FEmodelGeometry)
@ -156,6 +187,6 @@
endif
!$OMP END CRITICAL (write2out)
end subroutine
end subroutine FE_init
END MODULE FEsolving
end module FEsolving

28
code/config/debug.config
View File

@ -1,21 +1,19 @@
### $Id$ ###
### debugging parameters ###
verbosity 1 # level of detail of the debugging output (0-8)
# 0 : only version infos and all from "hypela2"/"umat"
# 1 : basic outputs from "CPFEM.f90", basic output from initialization routines, debug_info
# 2 : extensive outputs from "CPFEM.f90", extensive output from initialization routines
# 3 : basic outputs from "homogenization.f90"
# 4 : extensive outputs from "homogenization.f90"
# 5 : basic outputs from "crystallite.f90"
# 6 : extensive outputs from "crystallite.f90"
# 7 : basic outputs from the constitutive files
# 8 : extensive outputs from the constitutive files
selective 1 # >0 true to switch on e,i,g selective debugging
debug # debug.f90, possible keys: basic, extensive
math # math.f90, possible key: basic
FEsolving # FEsolving.f90, possible key: basic
math # math.f90, possible key: basic
material # material.f90, possible keys: basic, extensive
lattice # lattice.f90, possible key: basic
constitutive # constitutive_*.f90 possible keys: basic, extensive, selective
crystallite # crystallite.f90 possible keys: basic, extensive, selective
homogenization # homogenization_*.f90 possible keys: basic, extensive, selective
CPFEM # CPFEM.f90 possible keys: basic, selective
spectral # DAMASK_spectral.f90 possible keys: basic, fft, restart, divergence
#
# Parameters for selective
element 1 # selected element for debugging (synonymous: "el", "e")
ip 1 # selected integration point for debugging (synonymous: "integrationpoint", "i")
grain 1 # selected grain at ip for debugging (synonymous: "gr", "g")
### spectral solver debugging parameters ###
generalDebugSpectral 0 # > 0: general (algorithmical) debug outputs
divergenceDebugSpectral 0 # > 0: calculate more divergence measures and print them out

147
code/constitutive.f90
View File

@ -28,30 +28,36 @@
MODULE constitutive
!*** Include other modules ***
use prec
use prec, only: pInt, p_vec
implicit none
type(p_vec), dimension(:,:,:), allocatable :: &
constitutive_state0, & ! pointer array to microstructure at start of FE inc
constitutive_partionedState0, & ! pointer array to microstructure at start of homogenization inc
constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc
constitutive_state, & ! pointer array to current microstructure (end of converged time step)
constitutive_state_backup, & ! pointer array to backed up microstructure (end of converged time step)
constitutive_dotState, & ! pointer array to evolution of current microstructure
constitutive_previousDotState,& ! pointer array to previous evolution of current microstructure
constitutive_previousDotState2,& ! pointer array to 2nd previous evolution of current microstructure
constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure
constitutive_RK4dotState, & ! pointer array to evolution of microstructure defined by classical Runge-Kutta method
constitutive_aTolState ! pointer array to absolute state tolerance
type(p_vec), dimension(:,:,:), allocatable :: constitutive_state0, & ! pointer array to microstructure at start of FE inc
constitutive_partionedState0, & ! pointer array to microstructure at start of homogenization inc
constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc
constitutive_state, & ! pointer array to current microstructure (end of converged time step)
constitutive_state_backup, & ! pointer array to backed up microstructure (end of converged time step)
constitutive_dotState, & ! pointer array to evolution of current microstructure
constitutive_previousDotState,& ! pointer array to previous evolution of current microstructure
constitutive_previousDotState2,&! pointer array to 2nd previous evolution of current microstructure
constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure
constitutive_RK4dotState, & ! pointer array to evolution of microstructure defined by classical Runge-Kutta method
constitutive_aTolState ! pointer array to absolute state tolerance
type(p_vec), dimension(:,:,:,:), allocatable :: constitutive_RKCK45dotState ! pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method
integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array
constitutive_sizeState, & ! size of state array per grain
constitutive_sizePostResults ! size of postResults array per grain
integer(pInt) constitutive_maxSizeDotState, &
constitutive_maxSizeState, &
constitutive_maxSizePostResults
type(p_vec), dimension(:,:,:,:), allocatable :: &
constitutive_RKCK45dotState ! pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method
CONTAINS
integer(pInt), dimension(:,:,:), allocatable :: &
constitutive_sizeDotState, & ! size of dotState array
constitutive_sizeState, & ! size of state array per grain
constitutive_sizePostResults ! size of postResults array per grain
integer(pInt) :: &
constitutive_maxSizeDotState, &
constitutive_maxSizeState, &
constitutive_maxSizePostResults
contains
!****************************************
!* - constitutive_init
!* - constitutive_homogenizedC
@ -67,22 +73,36 @@ CONTAINS
!**************************************
!* Module initialization *
!**************************************
subroutine constitutive_init()
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt
use debug, only: debug_verbosity
use numerics, only: numerics_integrator
use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_titanmod
use constitutive_dislotwin
use constitutive_nonlocal
subroutine constitutive_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: debug_what, &
debug_constitutive, &
debug_levelBasic
use numerics, only: numerics_integrator
use IO, only: IO_error, &
IO_open_file, &
IO_open_jobFile_stat, &
IO_write_jobFile
use mesh, only: mesh_maxNips, &
mesh_NcpElems, &
mesh_element,FE_Nips
use material, only: material_phase, &
material_Nphase, &
material_localFileExt, &
material_configFile, &
phase_name, &
phase_constitution, &
phase_constitutionInstance, &
phase_Noutput, &
homogenization_Ngrains, &
homogenization_maxNgrains
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_titanmod
use constitutive_dislotwin
use constitutive_nonlocal
implicit none
integer(pInt), parameter :: fileunit = 200_pInt
integer(pInt) g, & ! grain number
i, & ! integration point number
@ -96,7 +116,7 @@ integer(pInt) g, & ! grain number
myNgrains
integer(pInt), dimension(:,:), pointer :: thisSize
character(len=64), dimension(:,:), pointer :: thisOutput
logical knownConstitution
logical :: knownConstitution
! --- PARSE CONSTITUTIONS FROM CONFIG FILE ---
@ -341,7 +361,7 @@ constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)
write(6,*) '<<<+- constitutive init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
if (debug_verbosity > 0_pInt) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt) then
write(6,'(a32,1x,7(i8,1x))') 'constitutive_state0: ', shape(constitutive_state0)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_subState0: ', shape(constitutive_subState0)
@ -371,17 +391,16 @@ function constitutive_homogenizedC(ipc,ip,el)
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt
use prec, only: pReal
use material, only: phase_constitution,material_phase
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_titanmod
use constitutive_dislotwin
use constitutive_nonlocal
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
implicit none
integer(pInt) :: ipc,ip,el
real(pReal), dimension(6,6) :: constitutive_homogenizedC
select case (phase_constitution(material_phase(ipc,ip,el)))
@ -415,17 +434,16 @@ function constitutive_averageBurgers(ipc,ip,el)
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt
use prec, only: pReal
use material, only: phase_constitution,material_phase
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_titanmod
use constitutive_dislotwin
use constitutive_nonlocal
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
implicit none
integer(pInt) :: ipc,ip,el
real(pReal) :: constitutive_averageBurgers
select case (phase_constitution(material_phase(ipc,ip,el)))
@ -456,7 +474,7 @@ endfunction
!* This function calculates from state needed variables *
!*********************************************************************
subroutine constitutive_microstructure(Temperature, Fe, Fp, ipc, ip, el)
use prec, only: pReal,pInt
use prec, only: pReal
use material, only: phase_constitution, &
material_phase
use constitutive_j2, only: constitutive_j2_label, &
@ -469,8 +487,8 @@ use constitutive_dislotwin, only: constitutive_dislotwin_label, &
constitutive_dislotwin_microstructure
use constitutive_nonlocal, only: constitutive_nonlocal_label, &
constitutive_nonlocal_microstructure
implicit none
implicit none
!*** input variables ***!
integer(pInt), intent(in):: ipc, & ! component-ID of current integration point
ip, & ! current integration point
@ -513,7 +531,7 @@ endsubroutine
!*********************************************************************
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ipc, ip, el)
use prec, only: pReal,pInt
use prec, only: pReal
use material, only: phase_constitution, &
material_phase
use constitutive_j2, only: constitutive_j2_label, &
@ -526,9 +544,8 @@ use constitutive_dislotwin, only: constitutive_dislotwin_label, &
constitutive_dislotwin_LpAndItsTangent
use constitutive_nonlocal, only: constitutive_nonlocal_label, &
constitutive_nonlocal_LpAndItsTangent
implicit none
!*** input variables ***!
integer(pInt), intent(in):: ipc, & ! component-ID of current integration point
ip, & ! current integration point
@ -573,10 +590,12 @@ endsubroutine
!*********************************************************************
subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, orientation, ipc, ip, el)
use prec, only: pReal, pInt
use prec, only: pReal, pLongInt
use debug, only: debug_cumDotStateCalls, &
debug_cumDotStateTicks, &
debug_verbosity
debug_what, &
debug_constitutive, &
debug_levelBasic
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: phase_constitution, &
@ -594,7 +613,6 @@ use constitutive_nonlocal, only: constitutive_nonlocal_dotState, &
constitutive_nonlocal_label
implicit none
!*** input variables
integer(pInt), intent(in) :: ipc, & ! component-ID of current integration point
ip, & ! current integration point
@ -608,15 +626,12 @@ real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),
orientation ! crystal orientation (quaternion)
real(pReal), dimension(6), intent(in) :: &
Tstar_v ! 2nd Piola Kirchhoff stress tensor (Mandel)
!*** output variables ***!
!*** local variables
integer(pLongInt) tick, tock, &
tickrate, &
maxticks
if (debug_verbosity > 0_pInt) then
if (iand(debug_what(debug_constitutive), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif
@ -640,7 +655,7 @@ select case (phase_constitution(material_phase(ipc,ip,el)))
end select
if (debug_verbosity > 6_pInt) then
if (iand(debug_what(debug_constitutive), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
!$OMP CRITICAL (debugTimingDotState)
debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt
@ -660,10 +675,12 @@ endsubroutine
!*********************************************************************
function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
use prec, only: pReal,pInt
use prec, only: pReal, pLongInt
use debug, only: debug_cumDotTemperatureCalls, &
debug_cumDotTemperatureTicks, &
debug_verbosity
debug_what, &
debug_constitutive, &
debug_levelBasic
use material, only: phase_constitution, &
material_phase
use constitutive_j2, only: constitutive_j2_dotTemperature, &
@ -676,8 +693,8 @@ use constitutive_dislotwin, only: constitutive_dislotwin_dotTemperature, &
constitutive_dislotwin_label
use constitutive_nonlocal, only: constitutive_nonlocal_dotTemperature, &
constitutive_nonlocal_label
implicit none
implicit none
!*** input variables
integer(pInt), intent(in) :: ipc, & ! component-ID of current integration point
ip, & ! current integration point
@ -695,7 +712,7 @@ integer(pLongInt) tick, tock, &
maxticks
if (debug_verbosity > 0_pInt) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt) then
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif
@ -718,7 +735,7 @@ select case (phase_constitution(material_phase(ipc,ip,el)))
end select
if (debug_verbosity > 6_pInt) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt) then
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
!$OMP CRITICAL (debugTimingDotTemperature)
debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt
@ -742,7 +759,7 @@ function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el)
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt
use prec, only: pReal
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: phase_constitution, &
@ -758,8 +775,8 @@ use constitutive_dislotwin, only: constitutive_dislotwin_postResults, &
constitutive_dislotwin_label
use constitutive_nonlocal, only: constitutive_nonlocal_postResults, &
constitutive_nonlocal_label
implicit none
implicit none
!*** input variables
integer(pInt), intent(in) :: ipc, & ! component-ID of current integration point
ip, & ! current integration point

242
code/constitutive_j2.f90
View File

@ -40,62 +40,79 @@
! tausat 63e6
! a 2.25
MODULE constitutive_j2
module constitutive_j2
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
private
character (len=*), parameter, public :: constitutive_j2_label = 'j2'
character (len=*), parameter :: constitutive_j2_label = 'j2'
integer(pInt), dimension(:), allocatable, public :: &
constitutive_j2_sizeDotState, &
constitutive_j2_sizeState, &
constitutive_j2_sizePostResults
integer(pInt), dimension(:,:), allocatable, target, public :: &
constitutive_j2_sizePostResult ! size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
constitutive_j2_output ! name of each post result output
integer(pInt), dimension(:), allocatable, private :: &
constitutive_j2_Noutput
real(pReal), dimension(:), allocatable, private ::&
constitutive_j2_C11, &
constitutive_j2_C12
real(pReal), dimension(:,:,:), allocatable, private :: &
constitutive_j2_Cslip_66
integer(pInt), dimension(:), allocatable :: constitutive_j2_sizeDotState, &
constitutive_j2_sizeState, &
constitutive_j2_sizePostResults
integer(pInt), dimension(:,:), allocatable,target :: constitutive_j2_sizePostResult ! size of each post result output
character(len=64), dimension(:,:), allocatable,target :: constitutive_j2_output ! name of each post result output
integer(pInt), dimension(:), allocatable :: constitutive_j2_Noutput
real(pReal), dimension(:), allocatable :: constitutive_j2_C11
real(pReal), dimension(:), allocatable :: constitutive_j2_C12
real(pReal), dimension(:,:,:), allocatable :: constitutive_j2_Cslip_66
!* Visco-plastic constitutive_j2 parameters
real(pReal), dimension(:), allocatable :: constitutive_j2_fTaylor
real(pReal), dimension(:), allocatable :: constitutive_j2_tau0
real(pReal), dimension(:), allocatable :: constitutive_j2_gdot0
real(pReal), dimension(:), allocatable :: constitutive_j2_n
real(pReal), dimension(:), allocatable :: constitutive_j2_h0
real(pReal), dimension(:), allocatable :: constitutive_j2_tausat
real(pReal), dimension(:), allocatable :: constitutive_j2_a
real(pReal), dimension(:), allocatable :: constitutive_j2_aTolResistance
real(pReal), dimension(:), allocatable, private :: &
constitutive_j2_fTaylor, &
constitutive_j2_tau0, &
constitutive_j2_gdot0, &
constitutive_j2_n, &
constitutive_j2_h0, &
constitutive_j2_tausat, &
constitutive_j2_a, &
constitutive_j2_aTolResistance
public :: constitutive_j2_init, &
constitutive_j2_stateInit, &
constitutive_j2_aTolState, &
constitutive_j2_homogenizedC, &
constitutive_j2_microstructure, &
constitutive_j2_LpAndItsTangent, &
constitutive_j2_dotState, &
constitutive_j2_dotTemperature, &
constitutive_j2_postResults
CONTAINS
!****************************************
!* - constitutive_j2_init
!* - constitutive_j2_stateInit
!* - constitutive_j2_homogenizedC
!* - constitutive_j2_microstructure
!* - constitutive_j2_LpAndItsTangent
!* - consistutive_j2_dotState
!* - consistutive_j2_postResults
!****************************************
contains
subroutine constitutive_j2_init(file)
subroutine constitutive_j2_init(myFile)
!**************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333
use IO
use material
use debug, only: debug_verbosity
integer(pInt), intent(in) :: file
use debug, only: debug_what, &
debug_constitutive, &
debug_levelBasic
implicit none
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 7_pInt
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k, mySize
character(len=64) tag
character(len=1024) line
integer(pInt) :: section = 0_pInt, maxNinstance, i,j,k, mySize
character(len=64) :: tag
character(len=1024) :: line
!$OMP CRITICAL (write2out)
write(6,*)
@ -107,41 +124,56 @@ subroutine constitutive_j2_init(file)
maxNinstance = int(count(phase_constitution == constitutive_j2_label),pInt)
if (maxNinstance == 0_pInt) return
if (debug_verbosity > 0_pInt) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a16,1x,i5)') '# instances:',maxNinstance
write(6,*)
!$OMP END CRITICAL (write2out)
endif
allocate(constitutive_j2_sizeDotState(maxNinstance)) ; constitutive_j2_sizeDotState = 0_pInt
allocate(constitutive_j2_sizeState(maxNinstance)) ; constitutive_j2_sizeState = 0_pInt
allocate(constitutive_j2_sizePostResults(maxNinstance)); constitutive_j2_sizePostResults = 0_pInt
allocate(constitutive_j2_sizePostResult(maxval(phase_Noutput), maxNinstance)); constitutive_j2_sizePostResult = 0_pInt
allocate(constitutive_j2_output(maxval(phase_Noutput), maxNinstance)) ; constitutive_j2_output = ''
allocate(constitutive_j2_Noutput(maxNinstance)) ; constitutive_j2_Noutput = 0_pInt
allocate(constitutive_j2_C11(maxNinstance)) ; constitutive_j2_C11 = 0.0_pReal
allocate(constitutive_j2_C12(maxNinstance)) ; constitutive_j2_C12 = 0.0_pReal
allocate(constitutive_j2_Cslip_66(6,6,maxNinstance)) ; constitutive_j2_Cslip_66 = 0.0_pReal
allocate(constitutive_j2_fTaylor(maxNinstance)) ; constitutive_j2_fTaylor = 0.0_pReal
allocate(constitutive_j2_tau0(maxNinstance)) ; constitutive_j2_tau0 = 0.0_pReal
allocate(constitutive_j2_gdot0(maxNinstance)) ; constitutive_j2_gdot0 = 0.0_pReal
allocate(constitutive_j2_n(maxNinstance)) ; constitutive_j2_n = 0.0_pReal
allocate(constitutive_j2_h0(maxNinstance)) ; constitutive_j2_h0 = 0.0_pReal
allocate(constitutive_j2_tausat(maxNinstance)) ; constitutive_j2_tausat = 0.0_pReal
allocate(constitutive_j2_a(maxNinstance)) ; constitutive_j2_a = 0.0_pReal
allocate(constitutive_j2_aTolResistance(maxNinstance)) ; constitutive_j2_aTolResistance = 0.0_pReal
allocate(constitutive_j2_sizeDotState(maxNinstance))
constitutive_j2_sizeDotState = 0_pInt
allocate(constitutive_j2_sizeState(maxNinstance))
constitutive_j2_sizeState = 0_pInt
allocate(constitutive_j2_sizePostResults(maxNinstance))
constitutive_j2_sizePostResults = 0_pInt
allocate(constitutive_j2_sizePostResult(maxval(phase_Noutput), maxNinstance))
constitutive_j2_sizePostResult = 0_pInt
allocate(constitutive_j2_output(maxval(phase_Noutput), maxNinstance))
constitutive_j2_output = ''
allocate(constitutive_j2_Noutput(maxNinstance))
constitutive_j2_Noutput = 0_pInt
allocate(constitutive_j2_C11(maxNinstance))
constitutive_j2_C11 = 0.0_pReal
allocate(constitutive_j2_C12(maxNinstance))
constitutive_j2_C12 = 0.0_pReal
allocate(constitutive_j2_Cslip_66(6,6,maxNinstance))
constitutive_j2_Cslip_66 = 0.0_pReal
allocate(constitutive_j2_fTaylor(maxNinstance))
constitutive_j2_fTaylor = 0.0_pReal
allocate(constitutive_j2_tau0(maxNinstance))
constitutive_j2_tau0 = 0.0_pReal
allocate(constitutive_j2_gdot0(maxNinstance))
constitutive_j2_gdot0 = 0.0_pReal
allocate(constitutive_j2_n(maxNinstance))
constitutive_j2_n = 0.0_pReal
allocate(constitutive_j2_h0(maxNinstance))
constitutive_j2_h0 = 0.0_pReal
allocate(constitutive_j2_tausat(maxNinstance))
constitutive_j2_tausat = 0.0_pReal
allocate(constitutive_j2_a(maxNinstance))
constitutive_j2_a = 0.0_pReal
allocate(constitutive_j2_aTolResistance(maxNinstance))
constitutive_j2_aTolResistance = 0.0_pReal
rewind(file)
line = ''
section = 0_pInt
rewind(myFile)
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
read(myFile,'(a1024)',END=100) line
enddo
do ! read thru sections of phase part
read(file,'(a1024)',END=100) line
read(myFile,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
@ -227,9 +259,7 @@ subroutine constitutive_j2_init(file)
enddo
return
endsubroutine
end subroutine constitutive_j2_init
!*********************************************************************
@ -237,16 +267,13 @@ endsubroutine
!*********************************************************************
pure function constitutive_j2_stateInit(myInstance)
use prec, only: pReal,pInt
implicit none
implicit none
integer(pInt), intent(in) :: myInstance
real(pReal), dimension(1) :: constitutive_j2_stateInit
integer(pInt), intent(in) :: myInstance
real(pReal), dimension(1) :: constitutive_j2_stateInit
constitutive_j2_stateInit = constitutive_j2_tau0(myInstance)
constitutive_j2_stateInit = constitutive_j2_tau0(myInstance)
return
endfunction
end function constitutive_j2_stateInit
!*********************************************************************
@ -254,22 +281,17 @@ endfunction
!*********************************************************************
pure function constitutive_j2_aTolState(myInstance)
use prec, only: pReal, &
pInt
implicit none
implicit none
!*** input variables
integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the constitution
!*** input variables
integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the constitution
!*** output variables
real(pReal), dimension(constitutive_j2_sizeState(myInstance)) :: &
!*** output variables
real(pReal), dimension(constitutive_j2_sizeState(myInstance)) :: &
constitutive_j2_aTolState ! relevant state values for the current instance of this constitution
!*** local variables
constitutive_j2_aTolState = constitutive_j2_aTolResistance(myInstance)
constitutive_j2_aTolState = constitutive_j2_aTolResistance(myInstance)
endfunction
end function constitutive_j2_aTolState
function constitutive_j2_homogenizedC(state,ipc,ip,el)
@ -281,22 +303,20 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el)
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use prec, only: p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
implicit none
integer(pInt), intent(in) :: ipc,ip,el
integer(pInt) matID
integer(pInt) :: matID
real(pReal), dimension(6,6) :: constitutive_j2_homogenizedC
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,matID)
return
endfunction
end function constitutive_j2_homogenizedC
subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el)
@ -308,11 +328,11 @@ subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el)
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use prec, only: p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el, matID
real(pReal) Temperature
@ -320,7 +340,7 @@ subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el)
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
endsubroutine
end subroutine constitutive_j2_microstructure
!****************************************************************
@ -329,9 +349,7 @@ endsubroutine
pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v, Temperature, state, g, ip, el)
!*** variables and functions from other modules ***!
use prec, only: pReal, &
pInt, &
p_vec
use prec, only: p_vec
use math, only: math_mul6x6, &
math_Mandel6to33, &
math_Plain3333to99
@ -342,7 +360,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v,
phase_constitutionInstance
implicit none
!*** input variables ***!
real(pReal), dimension(6), intent(in):: Tstar_dev_v ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in):: Temperature
@ -397,9 +414,7 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v,
dLp_dTstar_99 = math_Plain3333to99(gamma_dot / constitutive_j2_fTaylor(matID) * dLp_dTstar_3333 / norm_Tstar_dev)
end if
return
endsubroutine
end subroutine constitutive_j2_LpAndItsTangent
!****************************************************************
@ -408,9 +423,7 @@ endsubroutine
pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el)
!*** variables and functions from other modules ***!
use prec, only: pReal, &
pInt, &
p_vec
use prec, only: p_vec
use math, only: math_mul6x6
use mesh, only: mesh_NcpElems, &
mesh_maxNips
@ -419,7 +432,6 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el)
phase_constitutionInstance
implicit none
!*** input variables ***!
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: Temperature
@ -458,9 +470,7 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el)
! dotState
constitutive_j2_dotState = hardening * gamma_dot
return
endfunction
end function constitutive_j2_dotState
!****************************************************************
@ -469,11 +479,11 @@ endfunction
pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip, el)
!*** variables and functions from other modules ***!
use prec, only: pReal,pInt,p_vec
use prec, only: p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains
implicit none
implicit none
!*** input variables ***!
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: Temperature
@ -488,8 +498,7 @@ pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip,
! calculate dotTemperature
constitutive_j2_dotTemperature = 0.0_pReal
return
endfunction
end function constitutive_j2_dotTemperature
!*********************************************************************
@ -498,9 +507,7 @@ endfunction
pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip, el)
!*** variables and functions from other modules ***!
use prec, only: pReal, &
pInt, &
p_vec
use prec, only: p_vec
use math, only: math_mul6x6
use mesh, only: mesh_NcpElems, &
mesh_maxNips
@ -510,7 +517,6 @@ pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip
phase_Noutput
implicit none
!*** input variables ***!
real(pReal), dimension(6), intent(in):: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in):: Temperature, &
@ -561,8 +567,6 @@ pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip
end select
enddo
return
end function constitutive_j2_postResults
endfunction
END MODULE
end module constitutive_j2

63
code/constitutive_nonlocal.f90
View File

@ -155,7 +155,9 @@ use IO, only: IO_lc, &
IO_floatValue, &
IO_intValue, &
IO_error
use debug, only: debug_verbosity
use debug, only: debug_what, &
debug_constitutive, &
debug_levelBasic
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
FE_maxNipNeighbors
@ -212,7 +214,7 @@ character(len=1024) line
maxNinstance = int(count(phase_constitution == constitutive_nonlocal_label),pInt)
if (maxNinstance == 0) return ! we don't have to do anything if there's no instance for this constitutive law
if (debug_verbosity > 0) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a16,1x,i5)') '# instances:',maxNinstance
!$OMP END CRITICAL (write2out)
@ -894,8 +896,10 @@ use math, only: math_Mandel33to6, &
math_invert33, &
math_transpose33, &
pi
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_constitutive, &
debug_levelBasic, &
debug_levelSelective, &
debug_g, &
debug_i, &
debug_e
@ -1189,8 +1193,9 @@ state(g,ip,el)%p(12_pInt*ns+1:13_pInt*ns) = tauBack
#ifndef _OPENMP
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
.or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
.or. .not. iand(debug_what(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,g
write(6,*)
@ -1212,8 +1217,10 @@ subroutine constitutive_nonlocal_kinetics(v, tau, c, Temperature, state, g, ip,
use prec, only: pReal, &
pInt, &
p_vec
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_constitutive, &
debug_levelBasic, &
debug_levelSelective, &
debug_g, &
debug_i, &
debug_e
@ -1349,7 +1356,9 @@ endif
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
.or. .not. iand(debug_what(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_kinetics at el ip g',el,ip,g
write(6,*)
@ -1372,8 +1381,10 @@ use prec, only: pReal, &
p_vec
use math, only: math_Plain3333to99, &
math_mul6x6
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_constitutive, &
debug_levelBasic, &
debug_levelSelective, &
debug_g, &
debug_i, &
debug_e
@ -1491,8 +1502,9 @@ dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
#ifndef _OPENMP
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
.or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
.or. .not. iand(debug_what(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip g ',el,ip,g
write(6,*)
@ -1516,8 +1528,10 @@ use prec, only: pReal, &
DAMASK_NaN
use numerics, only: numerics_integrationMode
use IO, only: IO_error
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_constitutive, &
debug_levelBasic, &
debug_levelSelective, &
debug_g, &
debug_i, &
debug_e
@ -1628,8 +1642,9 @@ logical considerEnteringFlux, &
considerLeavingFlux
#ifndef _OPENMP
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) &
.or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
.or. .not. iand(debug_what(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g
write(6,*)
@ -1686,8 +1701,9 @@ forall (s = 1_pInt:ns, t = 1_pInt:4_pInt, rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pRea
gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgers(s,myInstance) * v(s,t)
#ifndef _OPENMP
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) &
.or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
.or. .not. iand(debug_what(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot
endif
@ -1700,7 +1716,7 @@ forall (s = 1_pInt:ns, t = 1_pInt:4_pInt, rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pRea
if (any(abs(gdot) > 0.0_pReal .and. 2.0_pReal * abs(v) * timestep > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then ! safety factor 2.0 (we use the reference volume and are for simplicity here)
#ifndef _OPENMP
if (debug_verbosity > 6_pInt) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt) then
write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ',maxval(abs(v)),' at a timestep of ',timestep
write(6,'(a)') '<< CONST >> enforcing cutback !!!'
@ -1966,7 +1982,7 @@ rhoDot = rhoDotFlux &
if ( any(rhoSgl(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < - constitutive_nonlocal_aTolRho(myInstance)) &
.or. any(rhoDip(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < - constitutive_nonlocal_aTolRho(myInstance))) then
#ifndef _OPENMP
if (debug_verbosity > 6_pInt) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt) then
write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
write(6,'(a)') '<< CONST >> enforcing cutback !!!'
endif
@ -1980,8 +1996,9 @@ endif
#ifndef _OPENMP
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
.or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
.or. .not. iand(debug_what(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation remobilization', rhoDotRemobilization(1:ns,1:8) * timestep
write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', rhoDotMultiplication(1:ns,1:4) * timestep
write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', rhoDotFlux(1:ns,1:8) * timestep

327
code/constitutive_phenopowerlaw.f90
View File

@ -70,92 +70,95 @@
!interaction_twintwin 1 1 1 1 1 1 1 1 10 10 10 10 10 10 10 10 10 10 10 10
!relevantResistance 1
MODULE constitutive_phenopowerlaw
module constitutive_phenopowerlaw
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
character (len=*), parameter :: &
constitutive_phenopowerlaw_label = 'phenopowerlaw'
character (len=*), parameter :: constitutive_phenopowerlaw_label = 'phenopowerlaw'
integer(pInt), dimension(:), allocatable :: &
constitutive_phenopowerlaw_sizeDotState, &
constitutive_phenopowerlaw_sizeState, &
constitutive_phenopowerlaw_sizePostResults, & ! cumulative size of post results
constitutive_phenopowerlaw_Noutput, & ! number of outputs per instance of this constitution
constitutive_phenopowerlaw_totalNslip, & ! no. of slip system used in simulation
constitutive_phenopowerlaw_totalNtwin, & ! no. of twin system used in simulation
constitutive_phenopowerlaw_structure
integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_sizeDotState, &
constitutive_phenopowerlaw_sizeState, &
constitutive_phenopowerlaw_sizePostResults ! cumulative size of post results
integer(pInt), dimension(:,:), allocatable,target :: constitutive_phenopowerlaw_sizePostResult ! size of each post result output
character(len=64), dimension(:,:), allocatable,target :: constitutive_phenopowerlaw_output ! name of each post result output
integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_Noutput ! number of outputs per instance of this constitution
integer(pInt), dimension(:,:), allocatable,target :: &
constitutive_phenopowerlaw_sizePostResult ! size of each post result output
character(len=32), dimension(:), allocatable :: constitutive_phenopowerlaw_structureName
integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_structure
integer(pInt), dimension(:,:), allocatable :: constitutive_phenopowerlaw_Nslip ! active number of slip systems per family
integer(pInt), dimension(:,:), allocatable :: constitutive_phenopowerlaw_Ntwin ! active number of twin systems per family
integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_totalNslip ! no. of slip system used in simulation
integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_totalNtwin ! no. of twin system used in simulation
integer(pInt), dimension(:,:), allocatable :: &
constitutive_phenopowerlaw_Nslip, & ! active number of slip systems per family
constitutive_phenopowerlaw_Ntwin ! active number of twin systems per family
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_CoverA
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C11
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C12
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C13
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C33
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C44
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_Cslip_66
!* Visco-plastic constitutive_phenomenological parameters
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_gdot0_slip
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_n_slip
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_tau0_slip
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_tausat_slip
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_gdot0_twin
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_n_twin
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_tau0_twin
character(len=64), dimension(:,:), allocatable,target :: &
constitutive_phenopowerlaw_output ! name of each post result output
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_spr
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinB
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinC
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinD
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinE
character(len=32), dimension(:), allocatable :: &
constitutive_phenopowerlaw_structureName
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_slipslip
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_sliptwin
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_twinslip
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_twintwin
real(pReal), dimension(:), allocatable :: &
constitutive_phenopowerlaw_CoverA, &
constitutive_phenopowerlaw_C11, &
constitutive_phenopowerlaw_C12, &
constitutive_phenopowerlaw_C13, &
constitutive_phenopowerlaw_C33, &
constitutive_phenopowerlaw_C44, &
constitutive_phenopowerlaw_gdot0_slip, &
constitutive_phenopowerlaw_n_slip, &
constitutive_phenopowerlaw_n_twin, &
constitutive_phenopowerlaw_gdot0_twin
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_slipslip
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_sliptwin
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_twinslip
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_twintwin
real(pReal), dimension(:,:), allocatable :: &
constitutive_phenopowerlaw_tau0_slip, &
constitutive_phenopowerlaw_tausat_slip, &
constitutive_phenopowerlaw_tau0_twin
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_slipslip
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_sliptwin
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_twinslip
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_twintwin
real(pReal), dimension(:), allocatable :: &
constitutive_phenopowerlaw_spr, &
constitutive_phenopowerlaw_twinB, &
constitutive_phenopowerlaw_twinC, &
constitutive_phenopowerlaw_twinD, &
constitutive_phenopowerlaw_twinE, &
constitutive_phenopowerlaw_h0_slipslip, &
constitutive_phenopowerlaw_h0_sliptwin, &
constitutive_phenopowerlaw_h0_twinslip, &
constitutive_phenopowerlaw_h0_twintwin, &
constitutive_phenopowerlaw_a_slip, &
constitutive_phenopowerlaw_aTolResistance
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_a_slip
real(pReal), dimension(:,:), allocatable :: &
constitutive_phenopowerlaw_interaction_slipslip, &
constitutive_phenopowerlaw_interaction_sliptwin, &
constitutive_phenopowerlaw_interaction_twinslip, &
constitutive_phenopowerlaw_interaction_twintwin
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_aTolResistance
real(pReal), dimension(:,:,:), allocatable :: &
constitutive_phenopowerlaw_hardeningMatrix_slipslip, &
constitutive_phenopowerlaw_hardeningMatrix_sliptwin, &
constitutive_phenopowerlaw_hardeningMatrix_twinslip, &
constitutive_phenopowerlaw_hardeningMatrix_twintwin, &
constitutive_phenopowerlaw_Cslip_66
CONTAINS
!****************************************
!* - constitutive_init
!* - constitutive_stateInit
!* - constitutive_homogenizedC
!* - constitutive_microstructure
!* - constitutive_LpAndItsTangent
!* - consistutive_dotState
!* - consistutive_postResults
!****************************************
public :: constitutive_phenopowerlaw_init
contains
subroutine constitutive_phenopowerlaw_init(myFile)
!**************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333
use math, only: math_Mandel3333to66, &
math_Voigt66to3333
use IO
use material
use debug, only: debug_verbosity
use debug, only: debug_what,&
debug_constitutive,&
debug_levelBasic
use lattice, only: lattice_initializeStructure, lattice_symmetryType, &
lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_maxNinteraction, lattice_NslipSystem, lattice_NtwinSystem, &
@ -164,13 +167,14 @@ subroutine constitutive_phenopowerlaw_init(myFile)
lattice_interactionTwinSlip, &
lattice_interactionTwinTwin
implicit none
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = lattice_maxNinteraction + 1_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k, f,o, &
mySize, myStructure, index_myFamily, index_otherFamily
character(len=64) tag
character(len=1024) line
character(len=64) :: tag
character(len=1024) :: line
!$OMP CRITICAL (write2out)
write(6,*)
@ -182,79 +186,96 @@ subroutine constitutive_phenopowerlaw_init(myFile)
maxNinstance = int(count(phase_constitution == constitutive_phenopowerlaw_label),pInt)
if (maxNinstance == 0) return
if (debug_verbosity > 0) then
if (iand(debug_what(debug_constitutive),debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a16,1x,i5)') '# instances:',maxNinstance
write(6,*)
!$OMP END CRITICAL (write2out)
endif
allocate(constitutive_phenopowerlaw_sizeDotState(maxNinstance)) ; constitutive_phenopowerlaw_sizeDotState = 0_pInt
allocate(constitutive_phenopowerlaw_sizeState(maxNinstance)) ; constitutive_phenopowerlaw_sizeState = 0_pInt
allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance)); constitutive_phenopowerlaw_sizePostResults = 0_pInt
allocate(constitutive_phenopowerlaw_sizePostResult(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_phenopowerlaw_sizePostResult = 0_pInt
allocate(constitutive_phenopowerlaw_output(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_phenopowerlaw_output = ''
allocate(constitutive_phenopowerlaw_Noutput(maxNinstance)) ; constitutive_phenopowerlaw_Noutput = 0_pInt
allocate(constitutive_phenopowerlaw_structureName(maxNinstance)) ; constitutive_phenopowerlaw_structureName = ''
allocate(constitutive_phenopowerlaw_structure(maxNinstance)) ; constitutive_phenopowerlaw_structure = 0_pInt
allocate(constitutive_phenopowerlaw_Nslip(lattice_maxNslipFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_Nslip = 0_pInt
allocate(constitutive_phenopowerlaw_Ntwin(lattice_maxNtwinFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_Ntwin = 0_pInt
allocate(constitutive_phenopowerlaw_totalNslip(maxNinstance)) ; constitutive_phenopowerlaw_totalNslip = 0_pInt !no. of slip system used in simulation (YJ.RO)
allocate(constitutive_phenopowerlaw_totalNtwin(maxNinstance)) ; constitutive_phenopowerlaw_totalNtwin = 0_pInt !no. of twin system used in simulation (YJ.RO)
allocate(constitutive_phenopowerlaw_CoverA(maxNinstance)) ; constitutive_phenopowerlaw_CoverA = 0.0_pReal
allocate(constitutive_phenopowerlaw_C11(maxNinstance)) ; constitutive_phenopowerlaw_C11 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C12(maxNinstance)) ; constitutive_phenopowerlaw_C12 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C13(maxNinstance)) ; constitutive_phenopowerlaw_C13 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C33(maxNinstance)) ; constitutive_phenopowerlaw_C33 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C44(maxNinstance)) ; constitutive_phenopowerlaw_C44 = 0.0_pReal
allocate(constitutive_phenopowerlaw_Cslip_66(6,6,maxNinstance)) ; constitutive_phenopowerlaw_Cslip_66 = 0.0_pReal
allocate(constitutive_phenopowerlaw_gdot0_slip(maxNinstance)) ; constitutive_phenopowerlaw_gdot0_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_n_slip(maxNinstance)) ; constitutive_phenopowerlaw_n_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_tau0_slip(lattice_maxNslipFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_tau0_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_tausat_slip(lattice_maxNslipFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_tausat_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_gdot0_twin(maxNinstance)) ; constitutive_phenopowerlaw_gdot0_twin = 0.0_pReal
allocate(constitutive_phenopowerlaw_n_twin(maxNinstance)) ; constitutive_phenopowerlaw_n_twin = 0.0_pReal
allocate(constitutive_phenopowerlaw_tau0_twin(lattice_maxNtwinFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_tau0_twin = 0.0_pReal
allocate(constitutive_phenopowerlaw_spr(maxNinstance)) ; constitutive_phenopowerlaw_spr = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinB(maxNinstance)) ; constitutive_phenopowerlaw_twinB = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinC(maxNinstance)) ; constitutive_phenopowerlaw_twinC = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinD(maxNinstance)) ; constitutive_phenopowerlaw_twinD = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinE(maxNinstance)) ; constitutive_phenopowerlaw_twinE = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_slipslip(maxNinstance)) ; constitutive_phenopowerlaw_h0_slipslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_sliptwin(maxNinstance)) ; constitutive_phenopowerlaw_h0_sliptwin = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_twinslip(maxNinstance)) ; constitutive_phenopowerlaw_h0_twinslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_twintwin(maxNinstance)) ; constitutive_phenopowerlaw_h0_twintwin = 0.0_pReal
allocate(constitutive_phenopowerlaw_sizeDotState(maxNinstance))
constitutive_phenopowerlaw_sizeDotState = 0_pInt
allocate(constitutive_phenopowerlaw_sizeState(maxNinstance))
constitutive_phenopowerlaw_sizeState = 0_pInt
allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance))
constitutive_phenopowerlaw_sizePostResults = 0_pInt
allocate(constitutive_phenopowerlaw_sizePostResult(maxval(phase_Noutput),maxNinstance))
constitutive_phenopowerlaw_sizePostResult = 0_pInt
allocate(constitutive_phenopowerlaw_output(maxval(phase_Noutput),maxNinstance))
constitutive_phenopowerlaw_output = ''
allocate(constitutive_phenopowerlaw_Noutput(maxNinstance))
constitutive_phenopowerlaw_Noutput = 0_pInt
allocate(constitutive_phenopowerlaw_structureName(maxNinstance))
constitutive_phenopowerlaw_structureName = ''
allocate(constitutive_phenopowerlaw_structure(maxNinstance))
constitutive_phenopowerlaw_structure = 0_pInt
allocate(constitutive_phenopowerlaw_Nslip(lattice_maxNslipFamily,maxNinstance))
constitutive_phenopowerlaw_Nslip = 0_pInt
allocate(constitutive_phenopowerlaw_Ntwin(lattice_maxNtwinFamily,maxNinstance))
constitutive_phenopowerlaw_Ntwin = 0_pInt
allocate(constitutive_phenopowerlaw_totalNslip(maxNinstance))
constitutive_phenopowerlaw_totalNslip = 0_pInt
allocate(constitutive_phenopowerlaw_totalNtwin(maxNinstance))
constitutive_phenopowerlaw_totalNtwin = 0_pInt
allocate(constitutive_phenopowerlaw_CoverA(maxNinstance))
constitutive_phenopowerlaw_CoverA = 0.0_pReal
allocate(constitutive_phenopowerlaw_C11(maxNinstance))
constitutive_phenopowerlaw_C11 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C12(maxNinstance))
constitutive_phenopowerlaw_C12 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C13(maxNinstance))
constitutive_phenopowerlaw_C13 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C33(maxNinstance))
constitutive_phenopowerlaw_C33 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C44(maxNinstance))
constitutive_phenopowerlaw_C44 = 0.0_pReal
allocate(constitutive_phenopowerlaw_Cslip_66(6,6,maxNinstance))
constitutive_phenopowerlaw_Cslip_66 = 0.0_pReal
allocate(constitutive_phenopowerlaw_gdot0_slip(maxNinstance))
constitutive_phenopowerlaw_gdot0_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_n_slip(maxNinstance))
constitutive_phenopowerlaw_n_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_tau0_slip(lattice_maxNslipFamily,maxNinstance))
constitutive_phenopowerlaw_tau0_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_tausat_slip(lattice_maxNslipFamily,maxNinstance))
constitutive_phenopowerlaw_tausat_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_gdot0_twin(maxNinstance))
constitutive_phenopowerlaw_gdot0_twin = 0.0_pReal
allocate(constitutive_phenopowerlaw_n_twin(maxNinstance))
constitutive_phenopowerlaw_n_twin = 0.0_pReal
allocate(constitutive_phenopowerlaw_tau0_twin(lattice_maxNtwinFamily,maxNinstance))
constitutive_phenopowerlaw_tau0_twin = 0.0_pReal
allocate(constitutive_phenopowerlaw_spr(maxNinstance))
constitutive_phenopowerlaw_spr = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinB(maxNinstance))
constitutive_phenopowerlaw_twinB = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinC(maxNinstance))
constitutive_phenopowerlaw_twinC = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinD(maxNinstance))
constitutive_phenopowerlaw_twinD = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinE(maxNinstance))
constitutive_phenopowerlaw_twinE = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_slipslip(maxNinstance))
constitutive_phenopowerlaw_h0_slipslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_sliptwin(maxNinstance))
constitutive_phenopowerlaw_h0_sliptwin = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_twinslip(maxNinstance))
constitutive_phenopowerlaw_h0_twinslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_twintwin(maxNinstance))
constitutive_phenopowerlaw_h0_twintwin = 0.0_pReal
allocate(constitutive_phenopowerlaw_interaction_slipslip(lattice_maxNinteraction,maxNinstance))
constitutive_phenopowerlaw_interaction_slipslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_interaction_sliptwin(lattice_maxNinteraction,maxNinstance))
constitutive_phenopowerlaw_interaction_sliptwin = 0.0_pReal
allocate(constitutive_phenopowerlaw_interaction_twinslip(lattice_maxNinteraction,maxNinstance))
constitutive_phenopowerlaw_interaction_twinslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_interaction_twintwin(lattice_maxNinteraction,maxNinstance))
constitutive_phenopowerlaw_interaction_slipslip = 0.0_pReal
constitutive_phenopowerlaw_interaction_sliptwin = 0.0_pReal
constitutive_phenopowerlaw_interaction_twinslip = 0.0_pReal
constitutive_phenopowerlaw_interaction_twintwin = 0.0_pReal
constitutive_phenopowerlaw_interaction_twintwin = 0.0_pReal
allocate(constitutive_phenopowerlaw_a_slip(maxNinstance))
constitutive_phenopowerlaw_a_slip = 0.0_pReal
constitutive_phenopowerlaw_a_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_aTolResistance(maxNinstance))
constitutive_phenopowerlaw_aTolResistance = 0.0_pReal
constitutive_phenopowerlaw_aTolResistance = 0.0_pReal
rewind(myFile)
line = ''
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
@ -525,20 +546,18 @@ subroutine constitutive_phenopowerlaw_init(myFile)
return
endsubroutine
end subroutine constitutive_phenopowerlaw_init
function constitutive_phenopowerlaw_stateInit(myInstance)
!*********************************************************************
!* initial microstructural state *
!*********************************************************************
use prec, only: pReal,pInt
use lattice, only: lattice_maxNslipFamily, lattice_maxNtwinFamily
implicit none
!* Definition of variables
implicit none
integer(pInt), intent(in) :: myInstance
integer(pInt) i
integer(pInt) :: i
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(myInstance)) :: constitutive_phenopowerlaw_stateInit
constitutive_phenopowerlaw_stateInit = 0.0_pReal
@ -559,7 +578,7 @@ function constitutive_phenopowerlaw_stateInit(myInstance)
enddo
return
endfunction
end function constitutive_phenopowerlaw_stateInit
!*********************************************************************
@ -567,10 +586,7 @@ endfunction
!*********************************************************************
pure function constitutive_phenopowerlaw_aTolState(myInstance)
use prec, only: pReal, &
pInt
implicit none
!*** input variables
integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the constitution
@ -582,7 +598,7 @@ real(pReal), dimension(constitutive_phenopowerlaw_sizeState(myInstance)) :: &
constitutive_phenopowerlaw_aTolState = constitutive_phenopowerlaw_aTolResistance(myInstance)
endfunction
end function constitutive_phenopowerlaw_aTolState
function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el)
@ -594,12 +610,11 @@ function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el)
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use prec, only: p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
implicit none
integer(pInt), intent(in) :: ipc,ip,el
integer(pInt) matID
real(pReal), dimension(6,6) :: constitutive_phenopowerlaw_homogenizedC
@ -610,7 +625,7 @@ function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el)
return
endfunction
end function constitutive_phenopowerlaw_homogenizedC
subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,ip,el)
@ -625,16 +640,15 @@ subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,ip,el
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
implicit none
integer(pInt) ipc,ip,el, matID
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
endsubroutine
end subroutine constitutive_phenopowerlaw_microstructure
subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el)
@ -649,7 +663,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
!* - Lp : plastic velocity gradient *
!* - dLp_dTstar : derivative of Lp (4th-rank tensor) *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use prec, only: p_vec
use math, only: math_Plain3333to99
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem
@ -657,8 +671,6 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
integer(pInt) matID,nSlip,nTwin,f,i,j,k,l,m,n, structID,index_Gamma,index_F,index_myFamily
real(pReal) Temperature
@ -741,7 +753,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
return
endsubroutine
end subroutine constitutive_phenopowerlaw_LpAndItsTangent
function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el)
@ -755,14 +767,13 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
!* OUTPUT: *
!* - constitutive_dotState : evolution of state variable *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use prec, only: p_vec
use lattice, only: lattice_Sslip_v, lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
implicit none
integer(pInt) ipc,ip,el
integer(pInt) matID,nSlip,nTwin,f,i,j, structID,index_Gamma,index_F,index_myFamily
real(pReal) Temperature,c_slipslip,c_sliptwin,c_twinslip,c_twintwin, ssat_offset
@ -864,9 +875,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
enddo
enddo
return
endfunction
end function constitutive_phenopowerlaw_dotState
!****************************************************************
@ -878,8 +887,8 @@ pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,stat
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems, mesh_maxNips
use material, only: homogenization_maxNgrains
implicit none
implicit none
!*** input variables ***!
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: Temperature
@ -894,8 +903,7 @@ pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,stat
! calculate dotTemperature
constitutive_phenopowerlaw_dotTemperature = 0.0_pReal
return
endfunction
end function constitutive_phenopowerlaw_dotTemperature
@ -914,9 +922,8 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
lattice_NslipSystem,lattice_NtwinSystem
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput
implicit none
!* Definition of variables
implicit none
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: dt,Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
@ -1006,8 +1013,6 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
end select
enddo
return
end function constitutive_phenopowerlaw_postResults
endfunction
END MODULE
end module constitutive_phenopowerlaw

328
code/constitutive_titanmod.f90
View File

@ -48,143 +48,175 @@ MODULE constitutive_titanmod
!* Include other modules
use prec, only: pReal,pInt
implicit none
!* Lists of states and physical parameters
character(len=*), parameter :: constitutive_titanmod_label = 'titanmod'
character(len=18), dimension(3), parameter:: constitutive_titanmod_listBasicSlipStates = (/'rho_edge ', &
'rho_screw ', &
'shear_system'/)
character(len=*), parameter :: &
constitutive_titanmod_label = 'titanmod'
character(len=18), dimension(3), parameter :: &
constitutive_titanmod_listBasicSlipStates = (/'rho_edge ', &
'rho_screw ', &
'shear_system'/)
character(len=18), dimension(1), parameter:: constitutive_titanmod_listBasicTwinStates = (/'gdot_twin'/)
character(len=18), dimension(1), parameter :: &
constitutive_titanmod_listBasicTwinStates = (/'gdot_twin'/)
character(len=19), dimension(11), parameter:: constitutive_titanmod_listDependentSlipStates =(/'segment_edge ', &
'segment_screw ', &
'resistance_edge ', &
'resistance_screw ', &
'tau_slip ', &
'velocity_edge ', &
'velocity_screw ', &
'gdot_slip_edge ', &
'gdot_slip_screw ', &
'stressratio_edge_p ', &
'stressratio_screw_p' &
/)
character(len=19), dimension(11), parameter :: &
constitutive_titanmod_listDependentSlipStates =(/'segment_edge ', &
'segment_screw ', &
'resistance_edge ', &
'resistance_screw ', &
'tau_slip ', &
'velocity_edge ', &
'velocity_screw ', &
'gdot_slip_edge ', &
'gdot_slip_screw ', &
'stressratio_edge_p ', &
'stressratio_screw_p' &
/)
character(len=18), dimension(2), parameter:: constitutive_titanmod_listDependentTwinStates =(/'twin_fraction', &
'tau_twin ' &
/)
character(len=18), dimension(2), parameter :: &
constitutive_titanmod_listDependentTwinStates =(/'twin_fraction', &
'tau_twin ' &
/)
real(pReal), parameter :: kB = 1.38e-23_pReal ! Boltzmann constant in J/Kelvin
!* Definition of global variables
integer(pInt), dimension(:), allocatable :: constitutive_titanmod_sizeDotState, & ! number of dotStates
constitutive_titanmod_sizeState, & ! total number of microstructural state variables
constitutive_titanmod_sizePostResults ! cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target :: constitutive_titanmod_sizePostResult ! size of each post result output
character(len=64), dimension(:,:), allocatable, target :: constitutive_titanmod_output ! name of each post result output
integer(pInt), dimension(:), allocatable :: constitutive_titanmod_Noutput ! number of outputs per instance of this constitution
character(len=32), dimension(:), allocatable :: constitutive_titanmod_structureName ! name of the lattice structure
integer(pInt), dimension(:), allocatable :: constitutive_titanmod_structure, & ! number representing the kind of lattice structure
constitutive_titanmod_totalNslip, & ! total number of active slip systems for each instance
constitutive_titanmod_totalNtwin ! total number of active twin systems for each instance
integer(pInt), dimension(:,:), allocatable :: constitutive_titanmod_Nslip, & ! number of active slip systems for each family and instance
constitutive_titanmod_Ntwin, & ! number of active twin systems for each family and instance
constitutive_titanmod_slipFamily, & ! lookup table relating active slip system to slip family for each instance
constitutive_titanmod_twinFamily, & ! lookup table relating active twin system to twin family for each instance
constitutive_titanmod_slipSystemLattice, & ! lookup table relating active slip system index to lattice slip system index for each instance
constitutive_titanmod_twinSystemLattice ! lookup table relating active twin system index to lattice twin system index for each instance
real(pReal), dimension(:), allocatable :: constitutive_titanmod_CoverA, & ! c/a ratio for hex type lattice
constitutive_titanmod_C11, & ! C11 element in elasticity matrix
constitutive_titanmod_C12, & ! C12 element in elasticity matrix
constitutive_titanmod_C13, & ! C13 element in elasticity matrix
constitutive_titanmod_C33, & ! C33 element in elasticity matrix
constitutive_titanmod_C44, & ! C44 element in elasticity matrix
constitutive_titanmod_debyefrequency, & !Debye frequency
constitutive_titanmod_kinkf0, & !Debye frequency
constitutive_titanmod_Gmod, & ! shear modulus
constitutive_titanmod_CAtomicVolume, & ! atomic volume in Bugers vector unit
constitutive_titanmod_dc, & ! prefactor for self-diffusion coefficient
constitutive_titanmod_twinhpconstant, & ! activation energy for dislocation climb
constitutive_titanmod_GrainSize, & ! grain size - Not being used
constitutive_titanmod_MaxTwinFraction, & ! maximum allowed total twin volume fraction
constitutive_titanmod_r, & ! r-exponent in twin nucleation rate
constitutive_titanmod_CEdgeDipMinDistance, & ! Not being used
constitutive_titanmod_Cmfptwin, & ! Not being used
constitutive_titanmod_Cthresholdtwin, & ! Not being used
constitutive_titanmod_aTolRho ! absolute tolerance for integration of dislocation density
real(pReal), dimension(:,:,:), allocatable :: constitutive_titanmod_Cslip_66 ! elasticity matrix in Mandel notation for each instance
real(pReal), dimension(:,:,:,:), allocatable :: constitutive_titanmod_Ctwin_66 ! twin elasticity matrix in Mandel notation for each instance
real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_titanmod_Cslip_3333 ! elasticity matrix for each instance
real(pReal), dimension(:,:,:,:,:,:), allocatable :: constitutive_titanmod_Ctwin_3333 ! twin elasticity matrix for each instance
real(pReal), dimension(:,:), allocatable :: constitutive_titanmod_rho_edge0, & ! initial edge dislocation density per slip system for each family and instance
constitutive_titanmod_rho_screw0, & ! initial screw dislocation density per slip system for each family and instance
constitutive_titanmod_shear_system0, & ! accumulated shear on each system
constitutive_titanmod_burgersPerSlipFamily, & ! absolute length of burgers vector [m] for each slip family and instance
constitutive_titanmod_burgersPerSlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance
constitutive_titanmod_burgersPerTwinFamily, & ! absolute length of burgers vector [m] for each twin family and instance
constitutive_titanmod_burgersPerTwinSystem, & ! absolute length of burgers vector [m] for each twin system and instance
constitutive_titanmod_f0_PerSlipFamily, & ! activation energy for glide [J] for each slip family and instance
constitutive_titanmod_f0_PerSlipSystem, & ! activation energy for glide [J] for each slip system and instance
constitutive_titanmod_twinf0_PerTwinFamily, & ! activation energy for glide [J] for each twin family and instance
constitutive_titanmod_twinf0_PerTwinSystem, & ! activation energy for glide [J] for each twin system and instance
constitutive_titanmod_twinshearconstant_PerTwinFamily, & ! activation energy for glide [J] for each twin family and instance
constitutive_titanmod_twinshearconstant_PerTwinSystem, & ! activation energy for glide [J] for each twin system and instance
constitutive_titanmod_tau0e_PerSlipFamily, & ! Initial yield stress for edge dislocations per slip family
constitutive_titanmod_tau0e_PerSlipSystem, & ! Initial yield stress for edge dislocations per slip system
constitutive_titanmod_tau0s_PerSlipFamily, & ! Initial yield stress for screw dislocations per slip family
constitutive_titanmod_tau0s_PerSlipSystem, & ! Initial yield stress for screw dislocations per slip system
constitutive_titanmod_twintau0_PerTwinFamily, & ! Initial yield stress for edge dislocations per twin family
constitutive_titanmod_twintau0_PerTwinSystem, & ! Initial yield stress for edge dislocations per twin system
constitutive_titanmod_capre_PerSlipFamily, & ! Capture radii for edge dislocations per slip family
constitutive_titanmod_capre_PerSlipSystem, & ! Capture radii for edge dislocations per slip system
constitutive_titanmod_caprs_PerSlipFamily, & ! Capture radii for screw dislocations per slip family
constitutive_titanmod_caprs_PerSlipSystem, & ! Capture radii for screw dislocations per slip system
constitutive_titanmod_pe_PerSlipFamily, & ! p-exponent in glide velocity
constitutive_titanmod_ps_PerSlipFamily, & ! p-exponent in glide velocity
constitutive_titanmod_qe_PerSlipFamily, & ! q-exponent in glide velocity
constitutive_titanmod_qs_PerSlipFamily, & ! q-exponent in glide velocity
constitutive_titanmod_pe_PerSlipSystem, & ! p-exponent in glide velocity
constitutive_titanmod_ps_PerSlipSystem, & ! p-exponent in glide velocity
constitutive_titanmod_qe_PerSlipSystem, & ! q-exponent in glide velocity
constitutive_titanmod_qs_PerSlipSystem, & ! q-exponent in glide velocity
constitutive_titanmod_twinp_PerTwinFamily, & ! p-exponent in glide velocity
constitutive_titanmod_twinq_PerTwinFamily, & ! q-exponent in glide velocity
constitutive_titanmod_twinp_PerTwinSystem, & ! p-exponent in glide velocity
constitutive_titanmod_twinq_PerTwinSystem, & ! p-exponent in glide velocity
constitutive_titanmod_v0e_PerSlipFamily, & ! edge dislocation velocity prefactor [m/s] for each family and instance
constitutive_titanmod_v0e_PerSlipSystem, & ! screw dislocation velocity prefactor [m/s] for each slip system and instance
constitutive_titanmod_v0s_PerSlipFamily, & ! edge dislocation velocity prefactor [m/s] for each family and instance
constitutive_titanmod_v0s_PerSlipSystem, & ! screw dislocation velocity prefactor [m/s] for each slip system and instance
constitutive_titanmod_twingamma0_PerTwinFamily, & ! edge dislocation velocity prefactor [m/s] for each family and instance
constitutive_titanmod_twingamma0_PerTwinSystem, & ! screw dislocation velocity prefactor [m/s] for each slip system and instance
constitutive_titanmod_kinkcriticallength_PerSlipFamily, & ! screw dislocation mobility prefactor for kink-pairs per slip family
constitutive_titanmod_kinkcriticallength_PerSlipSystem, & ! screw dislocation mobility prefactor for kink-pairs per slip system
constitutive_titanmod_twinsizePerTwinFamily, & ! twin thickness [m] for each twin family and instance
constitutive_titanmod_twinsizePerTwinSystem, & ! twin thickness [m] for each twin system and instance
constitutive_titanmod_CeLambdaSlipPerSlipFamily, & ! Adj. parameter for distance between 2 forest dislocations for each slip family and instance
constitutive_titanmod_CeLambdaSlipPerSlipSystem, & ! Adj. parameter for distance between 2 forest dislocations for each slip system and instance
constitutive_titanmod_CsLambdaSlipPerSlipFamily, & ! Adj. parameter for distance between 2 forest dislocations for each slip family and instance
constitutive_titanmod_CsLambdaSlipPerSlipSystem, & ! Adj. parameter for distance between 2 forest dislocations for each slip system and instance
constitutive_titanmod_twinLambdaSlipPerTwinFamily, & ! Adj. parameter for distance between 2 forest dislocations for each slip family and instance
constitutive_titanmod_twinLambdaSlipPerTwinSystem, & ! Adj. parameter for distance between 2 forest dislocations for each slip system and instance
constitutive_titanmod_interactionSlipSlip, & ! coefficients for slip-slip interaction for each interaction type and instance
constitutive_titanmod_interaction_ee, & ! coefficients for e-e interaction for each interaction type and instance
constitutive_titanmod_interaction_ss, & ! coefficients for s-s interaction for each interaction type and instance
constitutive_titanmod_interaction_es, & ! coefficients for e-s-twin interaction for each interaction type and instance
constitutive_titanmod_interactionSlipTwin, & ! coefficients for twin-slip interaction for each interaction type and instance
constitutive_titanmod_interactionTwinSlip, & ! coefficients for twin-slip interaction for each interaction type and instance
constitutive_titanmod_interactionTwinTwin ! coefficients for twin-twin interaction for each interaction type and instance
real(pReal), dimension(:,:,:), allocatable :: constitutive_titanmod_interactionMatrixSlipSlip, & ! interaction matrix of the different slip systems for each instance
constitutive_titanmod_interactionMatrix_ee, & ! interaction matrix of e-e for each instance
constitutive_titanmod_interactionMatrix_ss, & ! interaction matrix of s-s for each instance
constitutive_titanmod_interactionMatrix_es, & ! interaction matrix of e-s for each instance
constitutive_titanmod_interactionMatrixSlipTwin, & ! interaction matrix of slip systems with twin systems for each instance
constitutive_titanmod_interactionMatrixTwinSlip, & ! interaction matrix of twin systems with slip systems for each instance
constitutive_titanmod_interactionMatrixTwinTwin, & ! interaction matrix of the different twin systems for each instance
constitutive_titanmod_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance
constitutive_titanmod_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance
constitutive_titanmod_TwinforestProjectionEdge, & ! matrix of forest projections of edge dislocations in twin system for each instance
constitutive_titanmod_TwinforestProjectionScrew ! matrix of forest projections of screw dislocations in twin system for each instance
integer(pInt), dimension(:), allocatable :: &
constitutive_titanmod_sizeDotState, & ! number of dotStates
constitutive_titanmod_sizeState, & ! total number of microstructural state variables
constitutive_titanmod_sizePostResults ! cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target :: &
constitutive_titanmod_sizePostResult ! size of each post result output
character(len=64), dimension(:,:), allocatable, target :: &
constitutive_titanmod_output ! name of each post result output
integer(pInt), dimension(:), allocatable :: &
constitutive_titanmod_Noutput ! number of outputs per instance of this constitution
character(len=32), dimension(:), allocatable :: &
constitutive_titanmod_structureName ! name of the lattice structure
integer(pInt), dimension(:), allocatable :: &
constitutive_titanmod_structure, & ! number representing the kind of lattice structure
constitutive_titanmod_totalNslip, & ! total number of active slip systems for each instance
constitutive_titanmod_totalNtwin ! total number of active twin systems for each instance
integer(pInt), dimension(:,:), allocatable :: &
constitutive_titanmod_Nslip, & ! number of active slip systems for each family and instance
constitutive_titanmod_Ntwin, & ! number of active twin systems for each family and instance
constitutive_titanmod_slipFamily, & ! lookup table relating active slip system to slip family for each instance
constitutive_titanmod_twinFamily, & ! lookup table relating active twin system to twin family for each instance
constitutive_titanmod_slipSystemLattice, & ! lookup table relating active slip system index to lattice slip system index for each instance
constitutive_titanmod_twinSystemLattice ! lookup table relating active twin system index to lattice twin system index for each instance
real(pReal), dimension(:), allocatable :: &
constitutive_titanmod_CoverA, & ! c/a ratio for hex type lattice
constitutive_titanmod_C11, & ! C11 element in elasticity matrix
constitutive_titanmod_C12, & ! C12 element in elasticity matrix
constitutive_titanmod_C13, & ! C13 element in elasticity matrix
constitutive_titanmod_C33, & ! C33 element in elasticity matrix
constitutive_titanmod_C44, & ! C44 element in elasticity matrix
constitutive_titanmod_debyefrequency, & !Debye frequency
constitutive_titanmod_kinkf0, & !Debye frequency
constitutive_titanmod_Gmod, & ! shear modulus
constitutive_titanmod_CAtomicVolume, & ! atomic volume in Bugers vector unit
constitutive_titanmod_dc, & ! prefactor for self-diffusion coefficient
constitutive_titanmod_twinhpconstant, & ! activation energy for dislocation climb
constitutive_titanmod_GrainSize, & ! grain size - Not being used
constitutive_titanmod_MaxTwinFraction, & ! maximum allowed total twin volume fraction
constitutive_titanmod_r, & ! r-exponent in twin nucleation rate
constitutive_titanmod_CEdgeDipMinDistance, & ! Not being used
constitutive_titanmod_Cmfptwin, & ! Not being used
constitutive_titanmod_Cthresholdtwin, & ! Not being used
constitutive_titanmod_aTolRho ! absolute tolerance for integration of dislocation density
real(pReal), dimension(:,:,:), allocatable :: &
constitutive_titanmod_Cslip_66 ! elasticity matrix in Mandel notation for each instance
real(pReal), dimension(:,:,:,:), allocatable :: &
constitutive_titanmod_Ctwin_66 ! twin elasticity matrix in Mandel notation for each instance
real(pReal), dimension(:,:,:,:,:), allocatable :: &
constitutive_titanmod_Cslip_3333 ! elasticity matrix for each instance
real(pReal), dimension(:,:,:,:,:,:), allocatable :: &
constitutive_titanmod_Ctwin_3333 ! twin elasticity matrix for each instance
real(pReal), dimension(:,:), allocatable :: &
constitutive_titanmod_rho_edge0, & ! initial edge dislocation density per slip system for each family and instance
constitutive_titanmod_rho_screw0, & ! initial screw dislocation density per slip system for each family and instance
constitutive_titanmod_shear_system0, & ! accumulated shear on each system
constitutive_titanmod_burgersPerSlipFamily, & ! absolute length of burgers vector [m] for each slip family and instance
constitutive_titanmod_burgersPerSlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance
constitutive_titanmod_burgersPerTwinFamily, & ! absolute length of burgers vector [m] for each twin family and instance
constitutive_titanmod_burgersPerTwinSystem, & ! absolute length of burgers vector [m] for each twin system and instance
constitutive_titanmod_f0_PerSlipFamily, & ! activation energy for glide [J] for each slip family and instance
constitutive_titanmod_f0_PerSlipSystem, & ! activation energy for glide [J] for each slip system and instance
constitutive_titanmod_twinf0_PerTwinFamily, & ! activation energy for glide [J] for each twin family and instance
constitutive_titanmod_twinf0_PerTwinSystem, & ! activation energy for glide [J] for each twin system and instance
constitutive_titanmod_twinshearconstant_PerTwinFamily, & ! activation energy for glide [J] for each twin family and instance
constitutive_titanmod_twinshearconstant_PerTwinSystem, & ! activation energy for glide [J] for each twin system and instance
constitutive_titanmod_tau0e_PerSlipFamily, & ! Initial yield stress for edge dislocations per slip family
constitutive_titanmod_tau0e_PerSlipSystem, & ! Initial yield stress for edge dislocations per slip system
constitutive_titanmod_tau0s_PerSlipFamily, & ! Initial yield stress for screw dislocations per slip family
constitutive_titanmod_tau0s_PerSlipSystem, & ! Initial yield stress for screw dislocations per slip system
constitutive_titanmod_twintau0_PerTwinFamily, & ! Initial yield stress for edge dislocations per twin family
constitutive_titanmod_twintau0_PerTwinSystem, & ! Initial yield stress for edge dislocations per twin system
constitutive_titanmod_capre_PerSlipFamily, & ! Capture radii for edge dislocations per slip family
constitutive_titanmod_capre_PerSlipSystem, & ! Capture radii for edge dislocations per slip system
constitutive_titanmod_caprs_PerSlipFamily, & ! Capture radii for screw dislocations per slip family
constitutive_titanmod_caprs_PerSlipSystem, & ! Capture radii for screw dislocations per slip system
constitutive_titanmod_pe_PerSlipFamily, & ! p-exponent in glide velocity
constitutive_titanmod_ps_PerSlipFamily, & ! p-exponent in glide velocity
constitutive_titanmod_qe_PerSlipFamily, & ! q-exponent in glide velocity
constitutive_titanmod_qs_PerSlipFamily, & ! q-exponent in glide velocity
constitutive_titanmod_pe_PerSlipSystem, & ! p-exponent in glide velocity
constitutive_titanmod_ps_PerSlipSystem, & ! p-exponent in glide velocity
constitutive_titanmod_qe_PerSlipSystem, & ! q-exponent in glide velocity
constitutive_titanmod_qs_PerSlipSystem, & ! q-exponent in glide velocity
constitutive_titanmod_twinp_PerTwinFamily, & ! p-exponent in glide velocity
constitutive_titanmod_twinq_PerTwinFamily, & ! q-exponent in glide velocity
constitutive_titanmod_twinp_PerTwinSystem, & ! p-exponent in glide velocity
constitutive_titanmod_twinq_PerTwinSystem, & ! p-exponent in glide velocity
constitutive_titanmod_v0e_PerSlipFamily, & ! edge dislocation velocity prefactor [m/s] for each family and instance
constitutive_titanmod_v0e_PerSlipSystem, & ! screw dislocation velocity prefactor [m/s] for each slip system and instance
constitutive_titanmod_v0s_PerSlipFamily, & ! edge dislocation velocity prefactor [m/s] for each family and instance
constitutive_titanmod_v0s_PerSlipSystem, & ! screw dislocation velocity prefactor [m/s] for each slip system and instance
constitutive_titanmod_twingamma0_PerTwinFamily, & ! edge dislocation velocity prefactor [m/s] for each family and instance
constitutive_titanmod_twingamma0_PerTwinSystem, & ! screw dislocation velocity prefactor [m/s] for each slip system and instance
constitutive_titanmod_kinkcriticallength_PerSlipFamily, & ! screw dislocation mobility prefactor for kink-pairs per slip family
constitutive_titanmod_kinkcriticallength_PerSlipSystem, & ! screw dislocation mobility prefactor for kink-pairs per slip system
constitutive_titanmod_twinsizePerTwinFamily, & ! twin thickness [m] for each twin family and instance
constitutive_titanmod_twinsizePerTwinSystem, & ! twin thickness [m] for each twin system and instance
constitutive_titanmod_CeLambdaSlipPerSlipFamily, & ! Adj. parameter for distance between 2 forest dislocations for each slip family and instance
constitutive_titanmod_CeLambdaSlipPerSlipSystem, & ! Adj. parameter for distance between 2 forest dislocations for each slip system and instance
constitutive_titanmod_CsLambdaSlipPerSlipFamily, & ! Adj. parameter for distance between 2 forest dislocations for each slip family and instance
constitutive_titanmod_CsLambdaSlipPerSlipSystem, & ! Adj. parameter for distance between 2 forest dislocations for each slip system and instance
constitutive_titanmod_twinLambdaSlipPerTwinFamily, & ! Adj. parameter for distance between 2 forest dislocations for each slip family and instance
constitutive_titanmod_twinLambdaSlipPerTwinSystem, & ! Adj. parameter for distance between 2 forest dislocations for each slip system and instance
constitutive_titanmod_interactionSlipSlip, & ! coefficients for slip-slip interaction for each interaction type and instance
constitutive_titanmod_interaction_ee, & ! coefficients for e-e interaction for each interaction type and instance
constitutive_titanmod_interaction_ss, & ! coefficients for s-s interaction for each interaction type and instance
constitutive_titanmod_interaction_es, & ! coefficients for e-s-twin interaction for each interaction type and instance
constitutive_titanmod_interactionSlipTwin, & ! coefficients for twin-slip interaction for each interaction type and instance
constitutive_titanmod_interactionTwinSlip, & ! coefficients for twin-slip interaction for each interaction type and instance
constitutive_titanmod_interactionTwinTwin ! coefficients for twin-twin interaction for each interaction type and instance
real(pReal), dimension(:,:,:),allocatable :: &
constitutive_titanmod_interactionMatrixSlipSlip, & ! interaction matrix of the different slip systems for each instance
constitutive_titanmod_interactionMatrix_ee, & ! interaction matrix of e-e for each instance
constitutive_titanmod_interactionMatrix_ss, & ! interaction matrix of s-s for each instance
constitutive_titanmod_interactionMatrix_es, & ! interaction matrix of e-s for each instance
constitutive_titanmod_interactionMatrixSlipTwin, & ! interaction matrix of slip systems with twin systems for each instance
constitutive_titanmod_interactionMatrixTwinSlip, & ! interaction matrix of twin systems with slip systems for each instance
constitutive_titanmod_interactionMatrixTwinTwin, & ! interaction matrix of the different twin systems for each instance
constitutive_titanmod_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance
constitutive_titanmod_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance
constitutive_titanmod_TwinforestProjectionEdge, & ! matrix of forest projections of edge dislocations in twin system for each instance
constitutive_titanmod_TwinforestProjectionScrew ! matrix of forest projections of screw dislocations in twin system for each instance
CONTAINS
!****************************************
!* - constitutive_titanmod_init
@ -215,11 +247,11 @@ integer(pInt), intent(in) :: file
!* Local variables
integer(pInt), parameter :: maxNchunks = 21_pInt
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
integer(pInt) section,f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t,t1,t2,ns,nt,mySize,myStructure,maxTotalNslip, &
maxTotalNtwin
integer(pInt) :: section,f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t,t1,t2,ns,nt,&
mySize = 0_pInt,myStructure,maxTotalNslip,maxTotalNtwin
integer :: maxNinstance !no pInt
character(len=64) tag
character(len=1024) line
character(len=64) :: tag
character(len=1024) :: line
write(6,*)
write(6,*) '<<<+- constitutive_',trim(constitutive_titanmod_label),' init -+>>>'
@ -967,7 +999,7 @@ write(6,*) 'Determining elasticity matrix'
enddo
write(6,*) 'Init All done'
return
end subroutine
@ -977,8 +1009,8 @@ function constitutive_titanmod_stateInit(myInstance)
!*********************************************************************
use prec, only: pReal,pInt
use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily
implicit none
implicit none
!* Input-Output variables
integer(pInt) :: myInstance
real(pReal), dimension(constitutive_titanmod_sizeState(myInstance)) :: constitutive_titanmod_stateInit
@ -1062,7 +1094,6 @@ forall (t = 1_pInt:nt) &
resistance_twin0(t) = 0.0_pReal
constitutive_titanmod_stateInit(7_pInt*ns+nt+1_pInt:7_pInt*ns+2_pInt*nt)=resistance_twin0
return
end function
pure function constitutive_titanmod_aTolState(myInstance)
@ -1070,15 +1101,14 @@ pure function constitutive_titanmod_aTolState(myInstance)
!* absolute state tolerance *
!*********************************************************************
use prec, only: pReal, pInt
implicit none
implicit none
!* Input-Output variables
integer(pInt), intent(in) :: myInstance
real(pReal), dimension(constitutive_titanmod_sizeState(myInstance)) :: constitutive_titanmod_aTolState
constitutive_titanmod_aTolState = constitutive_titanmod_aTolRho(myInstance)
return
endfunction
pure function constitutive_titanmod_homogenizedC(state,g,ip,el)
@ -1092,8 +1122,8 @@ pure function constitutive_titanmod_homogenizedC(state,g,ip,el)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
implicit none
implicit none
!* Input-Output variables
integer(pInt), intent(in) :: g,ip,el
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
@ -1126,7 +1156,6 @@ do i=1_pInt,nt
enddo
return
end function
@ -1142,9 +1171,8 @@ subroutine constitutive_titanmod_microstructure(Temperature,state,g,ip,el)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
!use debug, only: debugger
implicit none
implicit none
!* Input-Output variables
integer(pInt), intent(in) :: g,ip,el
real(pReal), intent(in) :: Temperature
@ -1240,8 +1268,6 @@ forall (t = 1_pInt:nt) &
(dot_product((abs(state(g,ip,el)%p(2_pInt*ns+1_pInt:2_pInt*ns+nt))),&
constitutive_titanmod_interactionMatrixTwinTwin(1:nt,t,myInstance)))
return
end subroutine
@ -1265,8 +1291,8 @@ use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem, lattice_Stwin
implicit none
implicit none
!* Input-Output variables
integer(pInt), intent(in) :: g,ip,el
real(pReal), intent(in) :: Temperature
@ -1548,7 +1574,6 @@ dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
! write(6,'(a,/,9(9(f10.4,1x)/))') 'dLp_dTstar',dLp_dTstar
!endif
return
end subroutine
@ -1571,8 +1596,8 @@ use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem, lattice_Stwin_v
implicit none
implicit none
!* Input-Output variables
integer(pInt), intent(in) :: g,ip,el
real(pReal), intent(in) :: Temperature
@ -1696,8 +1721,6 @@ enddo
!write(6,'(a,/,4(3(f30.20,1x)/))') 'EdgeAnnihilation',DotRhoEdgeAnnihilation
!write(6,'(a,/,4(3(f30.20,1x)/))') 'ScrewAnnihilation',DotRhoScrewAnnihilation
return
end function
@ -1716,8 +1739,8 @@ pure function constitutive_titanmod_dotTemperature(Tstar_v,Temperature,state,g,i
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains
implicit none
implicit none
!* Input-Output variables
integer(pInt), intent(in) :: g,ip,el
real(pReal), intent(in) :: Temperature
@ -1727,7 +1750,6 @@ real(pReal) constitutive_titanmod_dotTemperature
constitutive_titanmod_dotTemperature = 0.0_pReal
return
end function
@ -1745,9 +1767,8 @@ pure function constitutive_titanmod_postResults(Tstar_v,Temperature,dt,state,g,i
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput
implicit none
!* Definition of variables
implicit none
integer(pInt), intent(in) :: g,ip,el
real(pReal), intent(in) :: dt,Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
@ -1898,7 +1919,6 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
end select
enddo
return
end function
END MODULE

446
code/crystallite.f90
View File

@ -33,8 +33,15 @@
MODULE crystallite
use prec, only: pReal, pInt
implicit none
implicit none
private :: crystallite_integrateStateFPI, &
crystallite_integrateStateEuler, &
crystallite_integrateStateAdaptiveEuler, &
crystallite_integrateStateRK4, &
crystallite_integrateStateRKCK45, &
crystallite_updateTemperature, &
crystallite_updateState
! ****************************************************************
! *** General variables for the crystallite calculation ***
! ****************************************************************
@ -104,11 +111,11 @@ subroutine crystallite_init(Temperature)
!*** variables and functions from other modules ***!
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt, &
pReal
use debug, only: debug_info, &
debug_reset, &
debug_verbosity
debug_what, &
debug_crystallite, &
debug_levelBasic
use math, only: math_I3, &
math_EulerToR, &
math_inv33, &
@ -383,7 +390,7 @@ call crystallite_stressAndItsTangent(.true.) ! request elastic
crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback
! *** Output to MARC output file ***
if (debug_verbosity > 0) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Temperature: ', shape(crystallite_Temperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dotTemperature: ', shape(crystallite_dotTemperature)
@ -435,10 +442,10 @@ if (debug_verbosity > 0) then
!$OMP END CRITICAL (write2out)
endif
call debug_info()
call debug_reset()
call debug_info
call debug_reset
endsubroutine
end subroutine crystallite_init
@ -448,102 +455,100 @@ endsubroutine
subroutine crystallite_stressAndItsTangent(updateJaco)
!*** variables and functions from other modules ***!
use prec, only: pInt, &
pReal
use numerics, only: subStepMinCryst, &
subStepSizeCryst, &
stepIncreaseCryst, &
pert_Fg, &
pert_method, &
nCryst, &
numerics_integrator, &
numerics_integrationMode, &
relevantStrain, &
Lp_frac, &
analyticJaco, &
time_sensitive
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_e, &
debug_i, &
debug_g, &
debug_CrystalliteLoopDistribution
use IO, only: IO_warning
use math, only: math_inv33, &
math_identity2nd, &
math_transpose33, &
math_mul33x33, &
math_mul66x6, &
math_Mandel6to33, &
math_Mandel33to6, &
math_I3, &
math_Plain3333to99, &
math_Plain99to3333, &
math_mul99x99, &
math_Mandel66to3333, &
math_mul3333xx33, &
math_invert, &
math_mul3333xx3333, &
math_spectralDecompositionSym33
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
use mesh, only: mesh_element, &
mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_Ngrains, &
homogenization_maxNgrains
use constitutive, only: constitutive_sizeState, &
constitutive_sizeDotState, &
constitutive_state, &
constitutive_state_backup, &
constitutive_subState0, &
constitutive_partionedState0, &
constitutive_homogenizedC, &
constitutive_dotState, &
constitutive_dotState_backup, &
constitutive_LpAndItsTangent
use numerics, only: subStepMinCryst, &
subStepSizeCryst, &
stepIncreaseCryst, &
pert_Fg, &
pert_method, &
nCryst, &
numerics_integrator, &
numerics_integrationMode, &
relevantStrain, &
Lp_frac, &
analyticJaco, &
time_sensitive
use debug, only: debug_what, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
debug_CrystalliteLoopDistribution
use IO, only: IO_warning
use math, only: math_inv33, &
math_identity2nd, &
math_transpose33, &
math_mul33x33, &
math_mul66x6, &
math_Mandel6to33, &
math_Mandel33to6, &
math_I3, &
math_Plain3333to99, &
math_Plain99to3333, &
math_mul99x99, &
math_Mandel66to3333, &
math_mul3333xx33, &
math_invert, &
math_mul3333xx3333, &
math_spectralDecompositionSym33
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
use mesh, only: mesh_element, &
mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_Ngrains, &
homogenization_maxNgrains
use constitutive, only: constitutive_sizeState, &
constitutive_sizeDotState, &
constitutive_state, &
constitutive_state_backup, &
constitutive_subState0, &
constitutive_partionedState0, &
constitutive_homogenizedC, &
constitutive_dotState, &
constitutive_dotState_backup, &
constitutive_LpAndItsTangent
implicit none
!*** input variables ***!
logical, intent(in) :: updateJaco ! flag indicating wehther we want to update the Jacobian (stiffness) or not
!*** output variables ***!
logical, intent(in) :: updateJaco ! flag indicating wehther we want to update the Jacobian (stiffness) or not
!*** local variables ***!
real(pReal) myPert, & ! perturbation with correct sign
formerSubStep
real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient
Fe_guess, & ! guess for elastic deformation gradient
Tstar ! 2nd Piola-Kirchhoff stress tensor
real(pReal) myPert, & ! perturbation with correct sign
formerSubStep
real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient
Fe_guess, & ! guess for elastic deformation gradient
Tstar ! 2nd Piola-Kirchhoff stress tensor
real(pReal), dimension(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
dPdF_perturbation1, &
dPdF_perturbation2
dPdF_perturbation1, &
dPdF_perturbation2
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
F_backup, &
Fp_backup, &
InvFp_backup, &
Fe_backup, &
Lp_backup, &
P_backup
F_backup, &
Fp_backup, &
InvFp_backup, &
Fe_backup, &
Lp_backup, &
P_backup
real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
Tstar_v_backup
Tstar_v_backup
real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
Temperature_backup
integer(pInt) NiterationCrystallite, & ! number of iterations in crystallite loop
e, & ! element index
i, & ! integration point index
g, & ! grain index
k, &
l, &
h, &
o, &
p, &
j, &
perturbation , & ! loop counter for forward,backward perturbation mode
myNgrains, &
mySizeState, &
mySizeDotState
Temperature_backup
integer(pInt) NiterationCrystallite, & ! number of iterations in crystallite loop
e, & ! element index
i, & ! integration point index
g, & ! grain index
k, &
l, &
h, &
o, &
p, &
j, &
perturbation , & ! loop counter for forward,backward perturbation mode
myNgrains, &
mySizeState, &
mySizeDotState
logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
convergenceFlag_backup
! local variables used for calculating analytic Jacobian
@ -590,7 +595,8 @@ logical :: error
! --+>> INITIALIZE TO STARTING CONDITION <<+--
if (debug_verbosity > 4 .and. debug_e > 0 .and. debug_e <= mesh_NcpElems &
if(iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt&
.and. debug_e > 0 .and. debug_e <= mesh_NcpElems &
.and. debug_i > 0 .and. debug_i <= mesh_maxNips &
.and. debug_g > 0 .and. debug_g <= homogenization_maxNgrains) then
!$OMP CRITICAL (write2out)
@ -651,8 +657,9 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
if (crystallite_converged(g,i,e)) then
#ifndef _OPENMP
if (debug_verbosity > 4 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite),debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,f12.8,a,f12.8,a)') '<< CRYST >> winding forward from ', &
crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', &
crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent'
@ -675,7 +682,7 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress
!$OMP FLUSH(crystallite_subF0)
elseif (formerSubStep > subStepMinCryst) then ! this crystallite just converged
if (debug_verbosity > 4_pInt) then
if (iand(debug_what(debug_crystallite),debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionCrystallite)
debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) = &
debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) + 1_pInt
@ -696,8 +703,9 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
! cant restore dotState here, since not yet calculated in first cutback after initialization
!$OMP FLUSH(crystallite_invFp)
#ifndef _OPENMP
if (debug_verbosity > 4_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite),debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,f12.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
crystallite_subStep(g,i,e)
write(6,*)
@ -761,8 +769,10 @@ enddo
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp)))
endif
#ifndef _OPENMP
if (debug_verbosity > 4 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if(iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite),debug_levelSelective) /= 0_pInt)) then
write (6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write (6,*)
write (6,'(a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal
@ -821,7 +831,7 @@ if(updateJaco) then
myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step
do k = 1,3; do l = 1,3 ! ...alter individual components
#ifndef _OPENMP
if (debug_verbosity> 5) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a,2(1x,i1),1x,a)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]'
write(6,*)
@ -1074,7 +1084,7 @@ if(updateJaco) then
endif
endif ! jacobian calculation
endsubroutine
end subroutine crystallite_stressAndItsTangent
@ -1088,8 +1098,11 @@ subroutine crystallite_integrateStateRK4(gg,ii,ee)
use prec, only: pInt, &
pReal
use numerics, only: numerics_integrationMode
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
@ -1248,8 +1261,10 @@ do n = 1_pInt,4_pInt
if (crystallite_integrateStress(g,i,e,timeStepFraction(n))) then ! fraction of original times step
if (n == 4) then ! final integration step
#ifndef _OPENMP
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
@ -1261,14 +1276,14 @@ do n = 1_pInt,4_pInt
#endif
crystallite_converged(g,i,e) = .true. ! ... converged per definition
crystallite_todo(g,i,e) = .false. ! ... integration done
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(n,numerics_integrationMode) = &
debug_StateLoopDistribution(n,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionState)
endif
endif
else ! broken stress integration
else ! broken stress integration
if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
@ -1329,7 +1344,7 @@ if (.not. singleRun) then
endif
endif
endsubroutine
end subroutine crystallite_integrateStateRK4
@ -1341,10 +1356,11 @@ endsubroutine
subroutine crystallite_integrateStateRKCK45(gg,ii,ee)
!*** variables and functions from other modules ***!
use prec, only: pInt, &
pReal
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
@ -1371,15 +1387,10 @@ use constitutive, only: constitutive_sizeDotState, &
constitutive_microstructure
implicit none
!*** input variables ***!
integer(pInt), optional, intent(in):: ee, & ! element index
ii, & ! integration point index
gg ! grain index
!*** output variables ***!
!*** local variables ***!
integer(pInt) e, & ! element index in element loop
i, & ! integration point index in ip loop
@ -1475,7 +1486,7 @@ endif
! --- FIRST RUNGE KUTTA STEP ---
#ifndef _OPENMP
if (debug_verbosity > 5) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,1x,i1)') '<< CRYST >> RUNGE KUTTA STEP',1
endif
#endif
@ -1611,8 +1622,8 @@ do n = 1_pInt,5_pInt
! --- dot state and RK dot state---
#ifndef _OPENMP
if (debug_verbosity > 5) then
write(6,'(a,1x,i1)') '<< CRYST >> RUNGE KUTTA STEP',n+1
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,1x,i1)') '<< CRYST >> RUNGE KUTTA STEP',n+1_pInt
endif
#endif
!$OMP DO
@ -1669,13 +1680,14 @@ relTemperatureResiduum = 0.0_pReal
! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION
! CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE
stateResiduum(1:mySizeDotState,g,i,e) = ( db(1) * constitutive_RKCK45dotState(1,g,i,e)%p(1:mySizeDotState) &
+ db(2) * constitutive_RKCK45dotState(2,g,i,e)%p(1:mySizeDotState) &
+ db(3) * constitutive_RKCK45dotState(3,g,i,e)%p(1:mySizeDotState) &
+ db(4) * constitutive_RKCK45dotState(4,g,i,e)%p(1:mySizeDotState) &
+ db(5) * constitutive_RKCK45dotState(5,g,i,e)%p(1:mySizeDotState) &
+ db(6) * constitutive_RKCK45dotState(6,g,i,e)%p(1:mySizeDotState)) &
* crystallite_subdt(g,i,e)
stateResiduum(1:mySizeDotState,g,i,e) = &
( db(1) * constitutive_RKCK45dotState(1,g,i,e)%p(1:mySizeDotState) &
+ db(2) * constitutive_RKCK45dotState(2,g,i,e)%p(1:mySizeDotState) &
+ db(3) * constitutive_RKCK45dotState(3,g,i,e)%p(1:mySizeDotState) &
+ db(4) * constitutive_RKCK45dotState(4,g,i,e)%p(1:mySizeDotState) &
+ db(5) * constitutive_RKCK45dotState(5,g,i,e)%p(1:mySizeDotState) &
+ db(6) * constitutive_RKCK45dotState(6,g,i,e)%p(1:mySizeDotState)) &
* crystallite_subdt(g,i,e)
temperatureResiduum(g,i,e) = ( db(1) * RKCK45dotTemperature(1,g,i,e) &
+ db(2) * RKCK45dotTemperature(2,g,i,e) &
+ db(3) * RKCK45dotTemperature(3,g,i,e) &
@ -1735,8 +1747,9 @@ relTemperatureResiduum = 0.0_pReal
.and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature )
#ifndef _OPENMP
if (debug_verbosity > 5_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt&
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i3,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> absolute residuum tolerance', &
@ -1776,7 +1789,7 @@ relTemperatureResiduum = 0.0_pReal
if (crystallite_integrateStress(g,i,e)) then
crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
crystallite_todo(g,i,e) = .false. ! ... integration done
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(6,numerics_integrationMode) =&
debug_StateLoopDistribution(6,numerics_integrationMode) + 1_pInt
@ -1798,7 +1811,7 @@ relTemperatureResiduum = 0.0_pReal
! --- nonlocal convergence check ---
#ifndef _OPENMP
if (debug_verbosity > 5) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged'
write(6,*)
endif
@ -1810,7 +1823,7 @@ if (.not. singleRun) then
endif
endsubroutine
end subroutine crystallite_integrateStateRKCK45
@ -1821,10 +1834,11 @@ endsubroutine
subroutine crystallite_integrateStateAdaptiveEuler(gg,ii,ee)
!*** variables and functions from other modules ***!
use prec, only: pInt, &
pReal
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
@ -1856,9 +1870,6 @@ implicit none
integer(pInt), optional, intent(in):: ee, & ! element index
ii, & ! integration point index
gg ! grain index
!*** output variables ***!
!*** local variables ***!
integer(pInt) e, & ! element index in element loop
i, & ! integration point index in ip loop
@ -2046,8 +2057,9 @@ relTemperatureResiduum = 0.0_pReal
!$OMP FLUSH(relStateResiduum,relTemperatureResiduum)
#ifndef _OPENMP
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> absolute residuum tolerance', &
@ -2071,7 +2083,7 @@ relTemperatureResiduum = 0.0_pReal
.and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature ) then
crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
crystallite_todo(g,i,e) = .false. ! ... integration done
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(2,numerics_integrationMode) = debug_StateLoopDistribution(2,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionState)
@ -2087,7 +2099,7 @@ relTemperatureResiduum = 0.0_pReal
! --- NONLOCAL CONVERGENCE CHECK ---
#ifndef _OPENMP
if (debug_verbosity > 5) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged'
write(6,*)
endif
@ -2098,7 +2110,7 @@ if (.not. singleRun) then
endif
endif
endsubroutine
end subroutine crystallite_integrateStateAdaptiveEuler
@ -2109,11 +2121,12 @@ endsubroutine
subroutine crystallite_integrateStateEuler(gg,ii,ee)
!*** variables and functions from other modules ***!
use prec, only: pInt, &
pReal
use numerics, only: numerics_integrationMode
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
@ -2132,14 +2145,10 @@ use constitutive, only: constitutive_sizeDotState, &
constitutive_microstructure
implicit none
!*** input variables ***!
integer(pInt), optional, intent(in):: ee, & ! element index
ii, & ! integration point index
gg ! grain index
!*** output variables ***!
!*** local variables ***!
integer(pInt) e, & ! element index in element loop
i, & ! integration point index in ip loop
@ -2220,8 +2229,10 @@ if (numerics_integrationMode < 2) then
crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
+ crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
#ifndef _OPENMP
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
@ -2256,7 +2267,7 @@ endif
if (crystallite_todo(g,i,e)) then
if (crystallite_integrateStress(g,i,e)) then
crystallite_converged(g,i,e) = .true.
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(1,numerics_integrationMode) = debug_StateLoopDistribution(1,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionState)
@ -2284,7 +2295,7 @@ if (.not. singleRun) then
endif
endif
endsubroutine
end subroutine crystallite_integrateStateEuler
@ -2296,9 +2307,10 @@ endsubroutine
subroutine crystallite_integrateStateFPI(gg,ii,ee)
!*** variables and functions from other modules ***!
use prec, only: pInt, &
pReal
use debug, only: debug_verbosity, &
use debug, only: debug_what,&
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_StateLoopDistribution
use numerics, only: nState, &
numerics_integrationMode
@ -2447,8 +2459,9 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
enddo; enddo; enddo
!$OMP ENDDO
#ifndef _OPENMP
if (debug_verbosity > 5) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration'
endif
#endif
@ -2500,7 +2513,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
elseif (stateConverged .and. temperatureConverged) then ! check (private) logicals "stateConverged" and "temperatureConverged" instead of (shared) "crystallite_converged", so no need to flush the "crystallite_converged" array
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(NiterationState,numerics_integrationMode) = &
debug_StateLoopDistribution(NiterationState,numerics_integrationMode) + 1_pInt
@ -2529,7 +2542,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP END PARALLEL
#ifndef _OPENMP
if (debug_verbosity > 5) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
' grains converged after state integration no. ', NiterationState
write(6,*)
@ -2547,7 +2560,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged
#ifndef _OPENMP
if (debug_verbosity > 5) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)),' grains converged after non-local check'
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after state integration no. ',&
NiterationState
@ -2557,7 +2570,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
enddo ! crystallite convergence loop
endsubroutine
end subroutine crystallite_integrateStateFPI
@ -2568,7 +2581,6 @@ endsubroutine
subroutine crystallite_updateState(done, converged, g, i, e)
!*** variables and functions from other modules ***!
use prec, only: pInt
use numerics, only: rTol_crystalliteState
use constitutive, only: constitutive_dotState, &
constitutive_previousDotState, &
@ -2577,8 +2589,11 @@ use constitutive, only: constitutive_dotState, &
constitutive_state, &
constitutive_aTolState, &
constitutive_microstructure
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g
@ -2618,7 +2633,7 @@ residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)
- dotState(1:mySize) * crystallite_subdt(g,i,e)
if (any(residuum /= residuum)) then ! if NaN occured then return without changing the state
#ifndef _OPENMP
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState encountered NaN at el ip g ',e,i,g
endif
#endif
@ -2634,7 +2649,9 @@ converged = all( abs(residuum) < constitutive_aTolState(g,i,e)%p(1:mySize) &
.or. abs(residuum) < rTol_crystalliteState * abs(state(1:mySize)) )
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
if (converged) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState converged at el ip g ',e,i,g
else
@ -2656,7 +2673,7 @@ endif
constitutive_dotState(g,i,e)%p(1:mySize) = dotState(1:mySize)
constitutive_state(g,i,e)%p(1:mySize) = state(1:mySize)
endsubroutine
end subroutine crystallite_updateState
@ -2664,13 +2681,14 @@ endsubroutine
! update the temperature of the grain
! and tell whether it has converged
!********************************************************************
subroutine crystallite_updateTemperature(done, converged, g, i, e)
subroutine crystallite_updateTemperature(done, converged, g, i, e)
!*** variables and functions from other modules ***!
use prec, only: pInt
use numerics, only: rTol_crystalliteTemperature
use constitutive, only: constitutive_dotTemperature
use debug, only: debug_verbosity
use debug, only: debug_what, &
debug_crystallite, &
debug_levelBasic
!*** input variables ***!
integer(pInt), intent(in):: e, & ! element index
i, & ! integration point index
@ -2698,7 +2716,7 @@ residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) &
* crystallite_subdt(g,i,e)
if (residuum /= residuum) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateTemperature encountered NaN at el ip g ',e,i,g
endif
#endif
@ -2714,7 +2732,7 @@ done = .true.
converged = ( crystallite_Temperature(g,i,e) == 0.0_pReal &
.or. abs(residuum) < rTol_crystalliteTemperature * crystallite_Temperature(g,i,e))
endsubroutine
end subroutine crystallite_updateTemperature
@ -2731,18 +2749,18 @@ function crystallite_integrateStress(&
)
!*** variables and functions from other modules ***!
use prec, only: pReal, &
pInt, &
pLongInt
use prec, only: pLongInt
use numerics, only: nStress, &
aTol_crystalliteStress, &
rTol_crystalliteStress, &
iJacoLpresiduum, &
relevantStrain, &
numerics_integrationMode
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_what, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
@ -2832,7 +2850,9 @@ integer(pLongInt) tick, &
crystallite_integrateStress = .false.
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress at el ip g ',e,i,g
endif
#endif
@ -2861,9 +2881,11 @@ Lpguess = crystallite_Lp(1:3,1:3,g,i,e) ! ... and tak
invFp_current = math_inv33(Fp_current)
if (all(invFp_current == 0.0_pReal)) then ! ... failed?
#ifndef _OPENMP
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on invFp_current inversion at el ip g ',e,i,g
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite), debug_levelSelective) > 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,*)
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new(1:3,1:3))
endif
@ -2896,7 +2918,7 @@ LpLoop: do
if (NiterationStress > nStress) then
#ifndef _OPENMP
if (debug_verbosity > 4_pInt) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress reached loop limit at el ip g ',e,i,g
write(6,*)
endif
@ -2919,11 +2941,11 @@ LpLoop: do
!* calculate plastic velocity gradient and its tangent according to constitutive law
if (debug_verbosity > 0) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e)
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
!$OMP CRITICAL (debugTimingLpTangent)
debug_cumLpCalls = debug_cumLpCalls + 1_pInt
@ -2934,7 +2956,9 @@ LpLoop: do
endif
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) &
if (iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt) &
.and. numerics_integrationMode == 1_pInt) then
write(6,'(a,i3)') '<< CRYST >> iteration ', NiterationStress
write(6,*)
@ -2961,7 +2985,7 @@ LpLoop: do
!* NaN occured at regular speed -> return
if (steplength >= steplength0 .and. any(residuum /= residuum)) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el ip g ',e,i,g,&
' ; iteration ', NiterationStress,&
' >> returning..!'
@ -3009,7 +3033,7 @@ LpLoop: do
!* something went wrong at accelerated speed? -> return to regular speed and try again
else
#ifndef _OPENMP
if (debug_verbosity > 5) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,i8,1x,i2,1x,i3,1x,a,i3)') '<< CRYST >> integrateStress encountered high-speed crash at el ip g ',e,i,g,&
'; iteration ', NiterationStress
endif
@ -3022,7 +3046,7 @@ LpLoop: do
steplength_max = steplength - 1.0_pReal ! limit acceleration
steplength = steplength0 ! grinding halt
jacoCounter = 0_pInt ! reset counter for Jacobian update (we want to do an update next time!)
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionLeapfrogBreak)
debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) = &
debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) + 1_pInt
@ -3046,10 +3070,11 @@ LpLoop: do
call math_invert(9_pInt,dR_dLp,inv_dR_dLp,dummy,error) ! invert dR/dLp --> dLp/dR
if (error) then
#ifndef _OPENMP
if (debug_verbosity > 4_pInt) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g
if (debug_verbosity > 5_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,*)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dR_dLp)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dLp',transpose(dT_dLp)
@ -3091,10 +3116,12 @@ invFp_new = invFp_new/math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize
call math_invert33(invFp_new,Fp_new,det,error)
if (error) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',&
e,i,g, ' ; iteration ', NiterationStress
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,*)
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new)
endif
@ -3124,7 +3151,9 @@ crystallite_invFp(1:3,1:3,g,i,e) = invFp_new
crystallite_integrateStress = .true.
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) &
if (iand(debug_what(debug_crystallite),debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_what(debug_crystallite), debug_levelSelective) /= 0_pInt) &
.and. numerics_integrationMode == 1_pInt) then
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Cauchy / MPa', &
@ -3135,25 +3164,24 @@ if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug
endif
#endif
if (debug_verbosity > 4) then
if (iand(debug_what(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) = &
debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1_pInt
!$OMP END CRITICAL (distributionStress)
endif
endfunction
end function crystallite_integrateStress
!********************************************************************
! calculates orientations and disorientations (in case of single grain ips)
!********************************************************************
subroutine crystallite_orientations()
subroutine crystallite_orientations
!*** variables and functions from other modules ***!
use prec, only: pInt, &
pReal
use math, only: math_pDecomposition, &
math_RtoQuaternion, &
math_QuaternionDisorientation, &
@ -3267,7 +3295,7 @@ logical error
enddo
!$OMP END PARALLEL DO
endsubroutine
end subroutine crystallite_orientations
@ -3282,8 +3310,6 @@ function crystallite_postResults(&
)
!*** variables and functions from other modules ***!
use prec, only: pInt, &
pReal
use math, only: math_QuaternionToEuler, &
math_QuaternionToAxisAngle, &
math_mul33x33, &
@ -3396,7 +3422,7 @@ function crystallite_postResults(&
dt, g, i, e)
c = c + constitutive_sizePostResults(g,i,e)
endfunction
end function crystallite_postResults
END MODULE

674
code/debug.f90
View File

@ -19,349 +19,423 @@
!##############################################################
!* $Id$
!##############################################################
MODULE debug
module debug
!##############################################################
use prec
use prec, only: pInt, pReal, pLongInt
implicit none
character(len=64), parameter :: debug_configFile = 'debug.config' ! name of configuration file
integer(pInt), parameter :: debug_spectralGeneral = 1_pInt, &
debug_spectralDivergence = 2_pInt, &
debug_spectralRestart = 4_pInt, &
debug_spectralFFTW = 8_pInt
implicit none
private
integer(pInt), dimension(:,:), allocatable :: debug_StressLoopDistribution
integer(pInt), dimension(:,:), allocatable :: debug_LeapfrogBreakDistribution
integer(pInt), dimension(:,:), allocatable :: debug_StateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution
integer(pLongInt) :: debug_cumLpTicks = 0_pLongInt
integer(pLongInt) :: debug_cumDotStateTicks = 0_pLongInt
integer(pLongInt) :: debug_cumDotTemperatureTicks = 0_pLongInt
integer(pInt) :: debug_cumLpCalls = 0_pInt
integer(pInt) :: debug_cumDotStateCalls = 0_pInt
integer(pInt) :: debug_cumDotTemperatureCalls = 0_pInt
integer(pInt) :: debug_e = 1_pInt
integer(pInt) :: debug_i = 1_pInt
integer(pInt) :: debug_g = 1_pInt
integer(pInt), dimension(2) :: debug_stressMaxLocation = 0_pInt
integer(pInt), dimension(2) :: debug_stressMinLocation = 0_pInt
integer(pInt), dimension(2) :: debug_jacobianMaxLocation = 0_pInt
integer(pInt), dimension(2) :: debug_jacobianMinLocation = 0_pInt
real(pReal) :: debug_stressMax
real(pReal) :: debug_stressMin
real(pReal) :: debug_jacobianMax
real(pReal) :: debug_jacobianMin
logical :: debug_selectiveDebugger = .true.
integer(pInt) :: debug_verbosity = 1_pInt
integer(pInt) :: debug_spectral = 0_pInt
integer(pInt), parameter, public :: &
debug_levelSelective = 2_pInt**0_pInt, &
debug_levelBasic = 2_pInt**1_pInt, &
debug_levelExtensive = 2_pInt**2_pInt
integer(pInt), parameter, private :: &
debug_maxForAll = debug_levelExtensive
integer(pInt), parameter, public :: &
debug_spectralRestart = debug_maxForAll*2_pInt**1_pInt, &
debug_spectralFFTW = debug_maxForAll*2_pInt**2_pInt, &
debug_spectralDivergence = debug_maxForAll*2_pInt**3_pInt
CONTAINS
integer(pInt), parameter, public :: &
debug_debug = 1_pInt, &
debug_math = 2_pInt, &
debug_FEsolving = 3_pInt, &
debug_mesh = 4_pInt, & ! stores debug level for mesh part of DAMASK
debug_material = 5_pInt, & ! stores debug level for material part of DAMASK
debug_lattice = 6_pInt, & ! stores debug level for lattice part of DAMASK
debug_constitutive = 7_pInt, & ! stores debug level for constitutive part of DAMASK
debug_crystallite = 8_pInt, &
debug_homogenization = 9_pInt, &
debug_CPFEM = 10_pInt, &
debug_spectral = 11_pInt
integer(pInt), dimension(11+2), public :: & ! 11 for specific, and 2 for "all" and "other"
debug_what = 0_pInt
integer(pInt), public :: &
debug_cumLpCalls = 0_pInt, &
debug_cumDotStateCalls = 0_pInt, &
debug_cumDotTemperatureCalls = 0_pInt, &
debug_e = 1_pInt, &
debug_i = 1_pInt, &
debug_g = 1_pInt
integer(pLongInt), public :: &
debug_cumLpTicks = 0_pLongInt, &
debug_cumDotStateTicks = 0_pLongInt, &
debug_cumDotTemperatureTicks = 0_pLongInt
integer(pInt), dimension(2), public :: &
debug_stressMaxLocation = 0_pInt, &
debug_stressMinLocation = 0_pInt, &
debug_jacobianMaxLocation = 0_pInt, &
debug_jacobianMinLocation = 0_pInt
integer(pInt), dimension(:), allocatable, public :: &
debug_CrystalliteLoopDistribution, &
debug_MaterialpointStateLoopDistribution, &
debug_MaterialpointLoopDistribution
integer(pInt), dimension(:,:), allocatable, public :: &
debug_StressLoopDistribution, &
debug_LeapfrogBreakDistribution, &
debug_StateLoopDistribution
real(pReal), public :: &
debug_stressMax = -huge(1.0_pReal), &
debug_stressMin = huge(1.0_pReal), &
debug_jacobianMax = -huge(1.0_pReal), &
debug_jacobianMin = huge(1.0_pReal)
character(len=64), parameter, private :: &
debug_configFile = 'debug.config' ! name of configuration file
public :: debug_init, &
debug_reset, &
debug_info
contains
!********************************************************************
! initialize the debugging capabilities
!********************************************************************
subroutine debug_init()
subroutine debug_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt
use numerics, only: nStress, &
nState, &
nCryst, &
nMPstate, &
nHomog
use IO, only: IO_error, &
IO_open_file_stat, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_lc, &
IO_floatValue, &
IO_intValue
implicit none
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use numerics, only: nStress, &
nState, &
nCryst, &
nMPstate, &
nHomog
use IO, only: IO_error, &
IO_open_file_stat, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_lc, &
IO_floatValue, &
IO_intValue
!*** input variables ***!
implicit none
integer(pInt), parameter :: fileunit = 300_pInt
integer(pInt), parameter :: maxNchunks = 6_pInt
!*** output variables ***!
integer(pInt) :: i, what
integer(pInt), dimension(1+2*maxNchunks) :: positions
character(len=64) :: tag
character(len=1024) :: line
!*** local variables ***!
integer(pInt), parameter :: fileunit = 300_pInt
integer(pInt), parameter :: maxNchunks = 2_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions
character(len=64) tag
character(len=1024) line
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- debug init -+>>>'
write(6,*) '$Id$'
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- debug init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
!$OMP END CRITICAL (write2out)
allocate(debug_StressLoopDistribution(nStress,2)) ; debug_StressLoopDistribution = 0_pInt
allocate(debug_LeapfrogBreakDistribution(nStress,2)) ; debug_LeapfrogBreakDistribution = 0_pInt
allocate(debug_StateLoopDistribution(nState,2)) ; debug_StateLoopDistribution = 0_pInt
allocate(debug_CrystalliteLoopDistribution(nCryst+1)) ; debug_CrystalliteLoopDistribution = 0_pInt
allocate(debug_MaterialpointStateLoopDistribution(nMPstate)) ; debug_MaterialpointStateLoopDistribution = 0_pInt
allocate(debug_MaterialpointLoopDistribution(nHomog+1)) ; debug_MaterialpointLoopDistribution = 0_pInt
allocate(debug_StressLoopDistribution(nStress,2))
debug_StressLoopDistribution = 0_pInt
allocate(debug_LeapfrogBreakDistribution(nStress,2))
debug_LeapfrogBreakDistribution = 0_pInt
allocate(debug_StateLoopDistribution(nState,2))
debug_StateLoopDistribution = 0_pInt
allocate(debug_CrystalliteLoopDistribution(nCryst+1))
debug_CrystalliteLoopDistribution = 0_pInt
allocate(debug_MaterialpointStateLoopDistribution(nMPstate))
debug_MaterialpointStateLoopDistribution = 0_pInt
allocate(debug_MaterialpointLoopDistribution(nHomog+1))
debug_MaterialpointLoopDistribution = 0_pInt
! try to open the config file
if(IO_open_file_stat(fileunit,debug_configFile)) then
line = ''
! read variables from config file and overwrite parameters
do
read(fileunit,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('element','e','el')
debug_e = IO_intValue(line,positions,2_pInt)
case ('integrationpoint','i','ip')
debug_i = IO_intValue(line,positions,2_pInt)
case ('grain','g','gr')
debug_g = IO_intValue(line,positions,2_pInt)
case ('selective')
debug_selectiveDebugger = IO_intValue(line,positions,2_pInt) > 0_pInt
case ('verbosity')
debug_verbosity = IO_intValue(line,positions,2_pInt)
case ('(spectral)')
select case(IO_lc(IO_stringValue(line,positions,2_pInt)))
case('general')
debug_spectral = ior(debug_spectral, debug_spectralGeneral)
case('divergence')
debug_spectral = ior(debug_spectral, debug_spectralDivergence)
case('restart')
debug_spectral = ior(debug_spectral, debug_spectralRestart)
case('fftw', 'fft')
debug_spectral = ior(debug_spectral, debug_spectralFFTW)
endselect
endselect
enddo
100 close(fileunit)
! try to open the config file
if(IO_open_file_stat(fileunit,debug_configFile)) then
if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,*) ' ... using values from config file'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
! read variables from config file and overwrite parameters
do
read(fileunit,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('element','e','el')
debug_e = IO_intValue(line,positions,2_pInt)
case ('integrationpoint','i','ip')
debug_i = IO_intValue(line,positions,2_pInt)
case ('grain','g','gr')
debug_g = IO_intValue(line,positions,2_pInt)
end select
! no config file, so we use standard values
else
what = 0_pInt
select case(tag)
case ('debug')
what = debug_debug
case ('math')
what = debug_math
case ('fesolving', 'fe')
what = debug_FEsolving
case ('mesh')
what = debug_mesh
case ('material')
what = debug_material
case ('lattice')
what = debug_lattice
case ('constitutive')
what = debug_constitutive
case ('crystallite')
what = debug_crystallite
case ('homogenization')
what = debug_homogenization
case ('cpfem')
what = debug_CPFEM
case ('spectral')
what = debug_spectral
case ('all')
what = 12_pInt
case ('other')
what = 13_pInt
end select
if(what /= 0) then
do i = 2_pInt, maxNchunks
select case(IO_lc(IO_stringValue(line,positions,i)))
case('basic')
debug_what(what) = ior(debug_what(what), debug_levelBasic)
case('extensive')
debug_what(what) = ior(debug_what(what), debug_levelExtensive)
case('selective')
debug_what(what) = ior(debug_what(what), debug_levelSelective)
case('restart')
debug_what(what) = ior(debug_what(what), debug_spectralRestart)
case('fft','fftw')
debug_what(what) = ior(debug_what(what), debug_spectralFFTW)
case('divergence')
debug_what(what) = ior(debug_what(what), debug_spectralDivergence)
end select
enddo
endif
enddo
100 close(fileunit)
if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,*) ' ... using standard values'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
do i = 1_pInt, 11_pInt
if(debug_what(i) == 0) debug_what(i) = ior(debug_what(i), debug_what(13))
debug_what(i) = ior(debug_what(i), debug_what(12))
enddo
endif
if (iand(debug_what(debug_debug),debug_levelBasic) /= 0) then
!$OMP CRITICAL (write2out)
write(6,*) 'using values from config file'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a24,1x,i1)') 'verbose: ',debug_verbosity
write(6,'(a24,1x,l1)') 'selective: ',debug_selectiveDebugger
!$OMP END CRITICAL (write2out)
endif
if (debug_selectiveDebugger) then
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a24,1x,i8)') 'element: ',debug_e
write(6,'(a24,1x,i8)') 'ip: ',debug_i
write(6,'(a24,1x,i8)') 'grain: ',debug_g
!$OMP END CRITICAL (write2out)
endif
else
debug_e = 0_pInt ! switch off selective debugging
debug_i = 0_pInt
debug_g = 0_pInt
endif
!$OMP CRITICAL (write2out) ! bitwise coded
if (iand(debug_spectral,debug_spectralGeneral) > 0_pInt) write(6,'(a)') ' spectral general debugging'
if (iand(debug_spectral,debug_spectralDivergence) > 0_pInt) write(6,'(a)') ' spectral divergence debugging'
if (iand(debug_spectral,debug_spectralRestart) > 0_pInt) write(6,'(a)') ' spectral restart debugging'
if (iand(debug_spectral,debug_spectralFFTW) > 0_pInt) write(6,'(a)') ' spectral FFTW debugging'
!$OMP END CRITICAL (write2out)
! no config file, so we use standard values
else
if (iand(debug_what(debug_debug),debug_levelBasic) /= 0) then
!$OMP CRITICAL (write2out)
write(6,*) 'using standard values'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
endif
endsubroutine
!output switched on (debug level for debug must be extensive)
if (iand(debug_what(debug_debug),debug_levelExtensive) /= 0) then
!$OMP CRITICAL (write2out)
do i = 1_pInt, 11_pInt
if(debug_what(i) /= 0) then
if(i == debug_debug) write(6,'(a)') 'Debug debugging:'
if(i == debug_math) write(6,'(a)') 'Math debugging:'
if(i == debug_FEsolving) write(6,'(a)') 'FEsolving debugging:'
if(i == debug_mesh) write(6,'(a)') 'Mesh debugging:'
if(i == debug_material) write(6,'(a)') 'Material debugging:'
if(i == debug_lattice) write(6,'(a)') 'Lattice debugging:'
if(i == debug_constitutive) write(6,'(a)') 'Constitutive debugging:'
if(i == debug_crystallite) write(6,'(a)') 'Crystallite debugging:'
if(i == debug_homogenization) write(6,'(a)') 'Homogenization debugging:'
if(i == debug_CPFEM) write(6,'(a)') 'CPFEM debugging:'
if(i == debug_spectral) write(6,'(a)') 'Spectral solver debugging:'
if(iand(debug_what(i),debug_levelBasic) /= 0) write(6,'(a)') ' basic'
if(iand(debug_what(i),debug_levelExtensive) /= 0) write(6,'(a)') ' extensive'
if(iand(debug_what(i),debug_levelSelective) /= 0) then
write(6,'(a)') 'selective on:'
write(6,'(a24,1x,i8)') 'element: ',debug_e
write(6,'(a24,1x,i8)') 'ip: ',debug_i
write(6,'(a24,1x,i8)') 'grain: ',debug_g
endif
if(iand(debug_what(i),debug_spectralRestart) /= 0) write(6,'(a)') ' restart'
if(iand(debug_what(i),debug_spectralFFTW) /= 0) write(6,'(a)') ' FFTW'
if(iand(debug_what(i),debug_spectralDivergence)/= 0) write(6,'(a)') ' divergence'
endif
enddo
!$OMP END CRITICAL (write2out)
endif
end subroutine debug_init
!********************************************************************
! reset debug distributions
!********************************************************************
subroutine debug_reset()
subroutine debug_reset
use prec
implicit none
implicit none
debug_StressLoopDistribution = 0_pInt ! initialize debugging data
debug_LeapfrogBreakDistribution = 0_pInt
debug_StateLoopDistribution = 0_pInt
debug_CrystalliteLoopDistribution = 0_pInt
debug_MaterialpointStateLoopDistribution = 0_pInt
debug_MaterialpointLoopDistribution = 0_pInt
debug_cumLpTicks = 0_pLongInt
debug_cumDotStateTicks = 0_pLongInt
debug_cumDotTemperatureTicks = 0_pLongInt
debug_cumLpCalls = 0_pInt
debug_cumDotStateCalls = 0_pInt
debug_cumDotTemperatureCalls = 0_pInt
debug_stressMaxLocation = 0_pInt
debug_stressMinLocation = 0_pInt
debug_jacobianMaxLocation = 0_pInt
debug_jacobianMinLocation = 0_pInt
debug_stressMax = -huge(1.0_pReal)
debug_stressMin = huge(1.0_pReal)
debug_jacobianMax = -huge(1.0_pReal)
debug_jacobianMin = huge(1.0_pReal)
debug_StressLoopDistribution = 0_pInt ! initialize debugging data
debug_LeapfrogBreakDistribution = 0_pInt
debug_StateLoopDistribution = 0_pInt
debug_CrystalliteLoopDistribution = 0_pInt
debug_MaterialpointStateLoopDistribution = 0_pInt
debug_MaterialpointLoopDistribution = 0_pInt
debug_cumLpTicks = 0_pLongInt
debug_cumDotStateTicks = 0_pLongInt
debug_cumDotTemperatureTicks = 0_pLongInt
debug_cumLpCalls = 0_pInt
debug_cumDotStateCalls = 0_pInt
debug_cumDotTemperatureCalls = 0_pInt
debug_stressMaxLocation = 0_pInt
debug_stressMinLocation = 0_pInt
debug_jacobianMaxLocation = 0_pInt
debug_jacobianMinLocation = 0_pInt
debug_stressMax = -huge(1.0_pReal)
debug_stressMin = huge(1.0_pReal)
debug_jacobianMax = -huge(1.0_pReal)
debug_jacobianMin = huge(1.0_pReal)
endsubroutine
end subroutine debug_reset
!********************************************************************
! write debug statements to standard out
!********************************************************************
subroutine debug_info()
subroutine debug_info
use prec
use numerics, only: nStress, &
nState, &
nCryst, &
nMPstate, &
nHomog
implicit none
use numerics, only: nStress, &
nState, &
nCryst, &
nMPstate, &
nHomog
integer(pInt) i,integral
integer(pLongInt) tickrate
implicit none
integer(pInt) :: i,integral
integer(pLongInt) :: tickrate
call system_clock(count_rate=tickrate)
call system_clock(count_rate=tickrate)
if (debug_verbosity > 4) then
!$OMP CRITICAL (write2out)
!$OMP CRITICAL (write2out)
if (iand(debug_what(debug_crystallite),debug_levelBasic) /= 0) then
write(6,*)
write(6,*) 'DEBUG Info (from previous cycle)'
write(6,*)
write(6,'(a33,1x,i12)') 'total calls to LpAndItsTangent :',debug_cumLpCalls
if (debug_cumLpCalls > 0_pInt) then
write(6,'(a33,1x,f12.3)') 'total CPU time/s :',real(debug_cumLpTicks,pReal)&
/real(tickrate,pReal)
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
real(debug_cumLpTicks,pReal)*1.0e6_pReal/real(tickrate,pReal)/real(debug_cumLpCalls,pReal)
endif
write(6,*)
write(6,'(a33,1x,i12)') 'total calls to collectDotState :',debug_cumDotStateCalls
if (debug_cumdotStateCalls > 0_pInt) then
write(6,'(a33,1x,f12.3)') 'total CPU time/s :',real(debug_cumDotStateTicks,pReal)&
/real(tickrate,pReal)
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
real(debug_cumDotStateTicks,pReal)*1.0e6_pReal/real(tickrate,pReal)&
/real(debug_cumDotStateCalls,pReal)
endif
write(6,*)
write(6,'(a33,1x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls
if (debug_cumdotTemperatureCalls > 0_pInt) then
write(6,'(a33,1x,f12.3)') 'total CPU time/s :',real(debug_cumDotTemperatureTicks,pReal)&
/real(tickrate,pReal)
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
real(debug_cumDotTemperatureTicks,pReal)*1.0e6_pReal/real(tickrate,pReal)&
/real(debug_cumDotTemperatureCalls,pReal)
endif
write(6,*)
write(6,*) 'DEBUG Info (from previous cycle)'
write(6,*)
write(6,'(a33,1x,i12)') 'total calls to LpAndItsTangent :',debug_cumLpCalls
if (debug_cumLpCalls > 0_pInt) then
write(6,'(a33,1x,f12.3)') 'total CPU time/s :',real(debug_cumLpTicks,pReal)&
/real(tickrate,pReal)
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
real(debug_cumLpTicks,pReal)*1.0e6_pReal/real(tickrate,pReal)/real(debug_cumLpCalls,pReal)
endif
write(6,*)
write(6,'(a33,1x,i12)') 'total calls to collectDotState :',debug_cumDotStateCalls
if (debug_cumdotStateCalls > 0_pInt) then
write(6,'(a33,1x,f12.3)') 'total CPU time/s :',real(debug_cumDotStateTicks,pReal)&
/real(tickrate,pReal)
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
real(debug_cumDotStateTicks,pReal)*1.0e6_pReal/real(tickrate,pReal)&
/real(debug_cumDotStateCalls,pReal)
endif
write(6,*)
write(6,'(a33,1x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls
if (debug_cumdotTemperatureCalls > 0_pInt) then
write(6,'(a33,1x,f12.3)') 'total CPU time/s :',real(debug_cumDotTemperatureTicks,pReal)&
/real(tickrate,pReal)
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
real(debug_cumDotTemperatureTicks,pReal)*1.0e6_pReal/real(tickrate,pReal)&
/real(debug_cumDotTemperatureCalls,pReal)
endif
integral = 0_pInt
write(6,*)
write(6,*)
write(6,*) 'distribution_StressLoop : stress frogbreak stiffness frogbreak'
do i=1_pInt,nStress
if (any(debug_StressLoopDistribution(i,:) /= 0_pInt ) .or. &
any(debug_LeapfrogBreakDistribution(i,:) /= 0_pInt ) ) then
integral = integral + i*debug_StressLoopDistribution(i,1) + i*debug_StressLoopDistribution(i,2)
write(6,'(i25,1x,i10,1x,i10,1x,i10,1x,i10)') i,debug_StressLoopDistribution(i,1),debug_LeapfrogBreakDistribution(i,1), &
debug_StressLoopDistribution(i,2),debug_LeapfrogBreakDistribution(i,2)
endif
enddo
write(6,'(a15,i10,1x,i10,12x,i10)') ' total',integral,&
sum(debug_StressLoopDistribution(:,1)), &
sum(debug_StressLoopDistribution(:,2))
integral = 0_pInt
write(6,*)
write(6,*)
write(6,*) 'distribution_StressLoop : stress frogbreak stiffness frogbreak'
do i=1_pInt,nStress
if (any(debug_StressLoopDistribution(i,:) /= 0_pInt ) .or. &
any(debug_LeapfrogBreakDistribution(i,:) /= 0_pInt ) ) then
integral = integral + i*debug_StressLoopDistribution(i,1) + i*debug_StressLoopDistribution(i,2)
write(6,'(i25,1x,i10,1x,i10,1x,i10,1x,i10)') i,debug_StressLoopDistribution(i,1),debug_LeapfrogBreakDistribution(i,1), &
debug_StressLoopDistribution(i,2),debug_LeapfrogBreakDistribution(i,2)
endif
enddo
write(6,'(a15,i10,1x,i10,12x,i10)') ' total',integral,&
sum(debug_StressLoopDistribution(:,1)), &
sum(debug_StressLoopDistribution(:,2))
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteStateLoop :'
do i=1_pInt,nState
if (any(debug_StateLoopDistribution(i,:) /= 0)) then
integral = integral + i*debug_StateLoopDistribution(i,1) + i*debug_StateLoopDistribution(i,2)
write(6,'(i25,1x,i10,12x,i10)') i,debug_StateLoopDistribution(i,1),debug_StateLoopDistribution(i,2)
endif
enddo
write(6,'(a15,i10,1x,i10,12x,i10)') ' total',integral,&
sum(debug_StateLoopDistribution(:,1)), &
sum(debug_StateLoopDistribution(:,2))
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteStateLoop :'
do i=1_pInt,nState
if (any(debug_StateLoopDistribution(i,:) /= 0)) then
integral = integral + i*debug_StateLoopDistribution(i,1) + i*debug_StateLoopDistribution(i,2)
write(6,'(i25,1x,i10,12x,i10)') i,debug_StateLoopDistribution(i,1),debug_StateLoopDistribution(i,2)
endif
enddo
write(6,'(a15,i10,1x,i10,12x,i10)') ' total',integral,&
sum(debug_StateLoopDistribution(:,1)), &
sum(debug_StateLoopDistribution(:,2))
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteCutbackLoop :'
do i=1_pInt,nCryst+1_pInt
if (debug_CrystalliteLoopDistribution(i) /= 0) then
integral = integral + i*debug_CrystalliteLoopDistribution(i)
if (i <= nCryst) then
write(6,'(i25,1x,i10)') i,debug_CrystalliteLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1_pInt,'+',debug_CrystalliteLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
endif
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteCutbackLoop :'
do i=1_pInt,nCryst+1_pInt
if (debug_CrystalliteLoopDistribution(i) /= 0) then
integral = integral + i*debug_CrystalliteLoopDistribution(i)
if (i <= nCryst) then
write(6,'(i25,1x,i10)') i,debug_CrystalliteLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1_pInt,'+',debug_CrystalliteLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
if (iand(debug_what(debug_homogenization),debug_levelBasic) /= 0) then
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointStateLoop :'
do i=1_pInt,nMPstate
if (debug_MaterialpointStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointStateLoopDistribution(i)
write(6,'(i25,1x,i10)') i,debug_MaterialpointStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution)
!$OMP END CRITICAL (write2out)
endif
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointCutbackLoop :'
do i=1_pInt,nHomog+1_pInt
if (debug_MaterialpointLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointLoopDistribution(i)
if (i <= nHomog) then
write(6,'(i25,1x,i10)') i,debug_MaterialpointLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1_pInt,'+',debug_MaterialpointLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
if (debug_verbosity > 2) then
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*)
write(6,*) 'Extreme values of returned stress and jacobian'
write(6,*)
write(6,'(a39)') ' value el ip'
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') 'stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') 'jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
write(6,*)
endif
!$OMP END CRITICAL (write2out)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointStateLoop :'
do i=1_pInt,nMPstate
if (debug_MaterialpointStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointStateLoopDistribution(i)
write(6,'(i25,1x,i10)') i,debug_MaterialpointStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution)
end subroutine debug_info
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointCutbackLoop :'
do i=1_pInt,nHomog+1_pInt
if (debug_MaterialpointLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointLoopDistribution(i)
if (i <= nHomog) then
write(6,'(i25,1x,i10)') i,debug_MaterialpointLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1_pInt,'+',debug_MaterialpointLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,*)
write(6,*)
write(6,*) 'Extreme values of returned stress and jacobian'
write(6,*)
write(6,'(a39)') ' value el ip'
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') 'stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') 'jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
write(6,*)
!$OMP END CRITICAL (write2out)
endif
endsubroutine
END MODULE debug
end module debug

34
code/homogenization.f90
View File

@ -73,9 +73,8 @@ CONTAINS
!**************************************
subroutine homogenization_init(Temperature)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pReal,pInt
use math, only: math_I3
use debug, only: debug_verbosity
use debug, only: debug_what, debug_homogenization, debug_levelBasic
use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material
@ -207,7 +206,7 @@ allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpEl
write(6,*) '<<<+- homogenization init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
if (debug_verbosity > 0) then
if (iand(debug_what(debug_homogenization), debug_levelBasic) /= 0_pInt) then
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state)
@ -249,8 +248,6 @@ subroutine materialpoint_stressAndItsTangent(&
dt & ! time increment
)
use prec, only: pInt, &
pReal
use numerics, only: subStepMinHomog, &
subStepSizeHomog, &
stepIncreaseHomog, &
@ -289,10 +286,12 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_converged, &
crystallite_stressAndItsTangent, &
crystallite_orientations
use debug, only: debug_verbosity, &
use debug, only: debug_what, &
debug_homogenization, &
debug_levelBasic, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_selectiveDebugger, &
debug_MaterialpointLoopDistribution, &
debug_MaterialpointStateLoopDistribution
use math, only: math_pDecomposition
@ -306,7 +305,8 @@ use debug, only: debug_verbosity, &
! ------ initialize to starting condition ------
if (debug_verbosity > 2 .and. debug_e > 0 .and. debug_e <= mesh_NcpElems .and. debug_i > 0 .and. debug_i <= mesh_maxNips) then
if (iand(debug_what(debug_homogenization), debug_levelBasic) /= 0_pInt &
.and. debug_e > 0 .and. debug_e <= mesh_NcpElems .and. debug_i > 0 .and. debug_i <= mesh_maxNips) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i
@ -358,7 +358,9 @@ use debug, only: debug_verbosity, &
if ( materialpoint_converged(i,e) ) then
#ifndef _OPENMP
if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_homogenization), debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_what(debug_homogenization),debug_levelSelective) /= 0_pInt)) then
write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,/)') '<< HOMOG >> winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent'
@ -388,7 +390,7 @@ use debug, only: debug_verbosity, &
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad
!$OMP FLUSH(materialpoint_subF0)
elseif (materialpoint_requested(i,e)) then ! this materialpoint just converged ! already at final time (??)
if (debug_verbosity > 2) then
if (iand(debug_what(debug_homogenization), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionHomog)
debug_MaterialpointLoopDistribution(min(nHomog+1,NiterationHomog)) = &
debug_MaterialpointLoopDistribution(min(nHomog+1,NiterationHomog)) + 1
@ -402,6 +404,7 @@ use debug, only: debug_verbosity, &
subStepSizeHomog * materialpoint_subStep(i,e) <= subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense and...
!$OMP CRITICAL (setTerminallyIll)
write(6,*) 'Integration point ', i,' at element ', e, ' terminally ill'
terminallyIll = .true. ! ...one kills all
!$OMP END CRITICAL (setTerminallyIll)
else ! cutback makes sense
@ -409,7 +412,9 @@ use debug, only: debug_verbosity, &
!$OMP FLUSH(materialpoint_subStep)
#ifndef _OPENMP
if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then
if (iand(debug_what(debug_homogenization), debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_what(debug_homogenization), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,1x,f12.8,/)') &
'<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',&
materialpoint_subStep(i,e)
@ -499,7 +504,7 @@ use debug, only: debug_verbosity, &
endif
!$OMP FLUSH(materialpoint_converged)
if (materialpoint_converged(i,e)) then
if (debug_verbosity > 2) then
if (iand(debug_what(debug_homogenization), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionMPState)
debug_MaterialpointStateLoopdistribution(NiterationMPstate) = &
debug_MaterialpointStateLoopdistribution(NiterationMPstate) + 1
@ -594,7 +599,6 @@ subroutine homogenization_partitionDeformation(&
el & ! element
)
use prec, only: pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_partionedF0,crystallite_partionedF
@ -635,7 +639,6 @@ function homogenization_updateState(&
ip, & ! integration point
el & ! element
)
use prec, only: pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_P,crystallite_dPdF,crystallite_partionedF,crystallite_partionedF0 ! modified <<<updated 31.07.2009>>>
@ -683,7 +686,6 @@ subroutine homogenization_averageStressAndItsTangent(&
ip, & ! integration point
el & ! element
)
use prec, only: pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_P,crystallite_dPdF
@ -725,7 +727,6 @@ subroutine homogenization_averageTemperature(&
ip, & ! integration point
el & ! element
)
use prec, only: pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_Temperature
@ -760,7 +761,6 @@ function homogenization_postResults(&
ip, & ! integration point
el & ! element
)
use prec, only: pReal,pInt
use mesh, only: mesh_element
use material, only: homogenization_type
use homogenization_isostrain

121
code/homogenization_RGC.f90
View File

@ -33,8 +33,8 @@ MODULE homogenization_RGC
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
implicit none
character (len=*), parameter :: homogenization_RGC_label = 'rgc'
integer(pInt), dimension(:), allocatable :: homogenization_RGC_sizeState, &
@ -67,12 +67,21 @@ subroutine homogenization_RGC_init(&
)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt, pReal
use debug, only: debug_verbosity
use math, only: math_Mandel3333to66, math_Voigt66to3333,math_I3,math_sampleRandomOri,math_EulerToR,inRad
use debug, only: debug_what, &
debug_homogenization, &
debug_levelBasic, &
debug_levelExtensive
use math, only: math_Mandel3333to66,&
math_Voigt66to3333, &
math_I3, &
math_sampleRandomOri,&
math_EulerToR,&
INRAD
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use IO
use material
implicit none
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 4_pInt
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
@ -170,7 +179,7 @@ subroutine homogenization_RGC_init(&
endif
enddo
100 if (debug_verbosity == 4_pInt) then
100 if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
do i = 1_pInt,maxNinstance
write(6,'(a15,1x,i4)') 'instance: ', i
@ -227,9 +236,8 @@ endsubroutine
!* initial homogenization state *
!*********************************************************************
function homogenization_RGC_stateInit(myInstance)
use prec, only: pReal,pInt
implicit none
implicit none
!* Definition of variables
integer(pInt), intent(in) :: myInstance
real(pReal), dimension(homogenization_RGC_sizeState(myInstance)) :: homogenization_RGC_stateInit
@ -253,8 +261,10 @@ subroutine homogenization_RGC_partitionDeformation(&
ip, & ! my integration point
el & ! my element
)
use prec, only: pReal,pInt,p_vec
use debug, only: debug_verbosity
use prec, only: p_vec
use debug, only: debug_what, &
debug_homogenization, &
debug_levelExtensive
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance
use FEsolving, only: theInc,cycleCounter
@ -277,7 +287,7 @@ subroutine homogenization_RGC_partitionDeformation(&
!* Debugging the overall deformation gradient
if (debug_verbosity == 4_pInt) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' =========='
write(6,'(1x,a32)')'Overall deformation gradient: '
@ -304,7 +314,7 @@ subroutine homogenization_RGC_partitionDeformation(&
F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! resulting relaxed deformation gradient
!* Debugging the grain deformation gradients
if (debug_verbosity == 4_pInt) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain
do i = 1_pInt,3_pInt
@ -338,7 +348,11 @@ function homogenization_RGC_updateState(&
)
use prec, only: pReal,pInt,p_vec
use debug, only: debug_verbosity, debug_e, debug_i
use debug, only: debug_what, &
debug_homogenization,&
debug_levelExtensive, &
debug_e, &
debug_i
use math, only: math_invert
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains,homogenization_typeInstance, &
@ -390,7 +404,7 @@ function homogenization_RGC_updateState(&
drelax = state%p(1:3_pInt*nIntFaceTot) - state0%p(1:3_pInt*nIntFaceTot)
!* Debugging the obtained state
if (debug_verbosity == 4_pInt) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Obtained state: '
do i = 1_pInt,3_pInt*nIntFaceTot
@ -407,7 +421,7 @@ function homogenization_RGC_updateState(&
call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID)
!* Debugging the mismatch, stress and penalties of grains
if (debug_verbosity == 4_pInt) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
do iGrain = 1_pInt,nGrain
write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e15.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain)
@ -456,7 +470,7 @@ function homogenization_RGC_updateState(&
enddo
!* Debugging the residual stress
if (debug_verbosity == 4_pInt) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum
write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1_pInt,3_pInt)
@ -474,7 +488,8 @@ function homogenization_RGC_updateState(&
residLoc = int(maxloc(abs(tract)),pInt) ! get the position of the maximum residual
!* Debugging the convergent criteria
if (debug_verbosity == 4_pInt .and. debug_e == el .and. debug_i == ip) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a)')' '
write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el
@ -491,7 +506,8 @@ function homogenization_RGC_updateState(&
if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then
homogenization_RGC_updateState = .true.
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a55)')'... done and happy'
write(6,*)' '
@ -521,7 +537,8 @@ function homogenization_RGC_updateState(&
state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) ! the average rate of relaxation vectors
state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors
if (debug_verbosity == 4_pInt .and. debug_e == el .and. debug_i == ip) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30,1x,e15.8)')'Constitutive work: ',constitutiveWork
write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/real(nGrain,pReal), &
@ -545,7 +562,8 @@ function homogenization_RGC_updateState(&
!* Try to restart when residual blows up exceeding maximum bound
homogenization_RGC_updateState = (/.true.,.false./) ! with direct cut-back
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a55)')'... broken'
write(6,*)' '
@ -559,7 +577,8 @@ function homogenization_RGC_updateState(&
!* Otherwise, proceed with computing the Jacobian and state update
else
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a55)')'... not yet done'
write(6,*)' '
@ -615,7 +634,7 @@ function homogenization_RGC_updateState(&
enddo
!* Debugging the global Jacobian matrix of stress tangent
if (debug_verbosity == 4_pInt) then
if (iand(debug_what(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of stress'
do i = 1_pInt,3_pInt*nIntFaceTot
@ -671,7 +690,7 @@ function homogenization_RGC_updateState(&
enddo
!* Debugging the global Jacobian matrix of penalty tangent
if (debug_verbosity == 4) then
if (iand(debug_what(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1_pInt,3_pInt*nIntFaceTot
@ -691,7 +710,7 @@ function homogenization_RGC_updateState(&
! only in the main diagonal term
!* Debugging the global Jacobian matrix of numerical viscosity tangent
if (debug_verbosity == 4_pInt) then
if (iand(debug_what(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1_pInt,3_pInt*nIntFaceTot
@ -705,7 +724,7 @@ function homogenization_RGC_updateState(&
!* The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix
if (debug_verbosity == 4) then
if (iand(debug_what(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix (total)'
do i = 1_pInt,3_pInt*nIntFaceTot
@ -724,7 +743,7 @@ function homogenization_RGC_updateState(&
call math_invert(3_pInt*nIntFaceTot,jmatrix,jnverse,ival,error) ! Compute the inverse of the overall Jacobian matrix
!* Debugging the inverse Jacobian matrix
if (debug_verbosity == 4_pInt) then
if (iand(debug_what(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian inverse'
do i = 1_pInt,3_pInt*nIntFaceTot
@ -754,7 +773,7 @@ function homogenization_RGC_updateState(&
endif
!* Debugging the return state
if (debug_verbosity == 4_pInt) then
if (iand(debug_homogenization, debug_levelExtensive) > 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Returned state: '
do i = 1_pInt,3_pInt*nIntFaceTot
@ -784,13 +803,14 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
)
use prec, only: pReal,pInt,p_vec
use debug, only: debug_verbosity
use debug, only: debug_what, &
debug_homogenization,&
debug_levelExtensive
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance
use math, only: math_Plain3333to99
implicit none
!* Definition of variables
implicit none
real(pReal), dimension (3,3), intent(out) :: avgP
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P
@ -804,7 +824,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
Ngrains = homogenization_Ngrains(mesh_element(3,el))
!* Debugging the grain tangent
if (debug_verbosity == 4_pInt) then
if (iand(debug_what(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
do iGrain = 1_pInt,Ngrains
dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain))
@ -836,9 +856,8 @@ function homogenization_RGC_averageTemperature(&
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains, homogenization_Ngrains
implicit none
!* Definition of variables
implicit none
real(pReal), dimension (homogenization_maxNgrains), intent(in) :: Temperature
integer(pInt), intent(in) :: ip,el
real(pReal) homogenization_RGC_averageTemperature
@ -862,9 +881,8 @@ pure function homogenization_RGC_postResults(&
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use material, only: homogenization_typeInstance,homogenization_Noutput
implicit none
!* Definition of variables
implicit none
type(p_vec), intent(in) :: state
integer(pInt), intent(in) :: ip,el
!
@ -925,9 +943,8 @@ subroutine homogenization_RGC_stressPenalty(&
use math, only: math_civita,math_invert33
use material, only: homogenization_maxNgrains,homogenization_Ngrains
use numerics, only: xSmoo_RGC
implicit none
!* Definition of variables
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen
real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef
@ -1059,8 +1076,6 @@ subroutine homogenization_RGC_volumePenalty(&
use numerics, only: maxVolDiscr_RGC,volDiscrMod_RGC,volDiscrPow_RGC
implicit none
!* Definition of variables
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: vPen
real(pReal), intent(out) :: vDiscrep
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef
@ -1109,9 +1124,8 @@ function homogenization_RGC_surfaceCorrection(&
use prec, only: pReal,pInt,p_vec
use math, only: math_invert33,math_mul33x33
implicit none
!* Definition of variables
implicit none
real(pReal), dimension(3,3), intent(in) :: avgF
real(pReal), dimension(3) :: homogenization_RGC_surfaceCorrection
integer(pInt), intent(in) :: ip,el
@ -1154,9 +1168,8 @@ function homogenization_RGC_equivalentModuli(&
use prec, only: pReal,pInt,p_vec
use constitutive, only: constitutive_homogenizedC,constitutive_averageBurgers
implicit none
!* Definition of variables
implicit none
integer(pInt), intent(in) :: grainID,ip,el
real(pReal), dimension (6,6) :: elasTens
real(pReal), dimension(2) :: homogenization_RGC_equivalentModuli
@ -1186,9 +1199,8 @@ function homogenization_RGC_relaxationVector(&
)
use prec, only: pReal,pInt,p_vec
implicit none
!* Definition of variables
implicit none
real(pReal), dimension (3) :: homogenization_RGC_relaxationVector
integer(pInt), dimension (4), intent(in) :: intFace
type(p_vec), intent(in) :: state
@ -1215,9 +1227,8 @@ function homogenization_RGC_interfaceNormal(&
use prec, only: pReal,pInt,p_vec
use math, only: math_mul33x3
implicit none
!* Definition of variables
implicit none
real(pReal), dimension (3) :: homogenization_RGC_interfaceNormal
integer(pInt), dimension (4), intent(in) :: intFace
integer(pInt), intent(in) :: ip,el
@ -1249,9 +1260,8 @@ function homogenization_RGC_getInterface(&
iGrain3 & ! grain ID in 3D array
)
use prec, only: pReal,pInt,p_vec
implicit none
!* Definition of variables
implicit none
integer(pInt), dimension (4) :: homogenization_RGC_getInterface
integer(pInt), dimension (3), intent(in) :: iGrain3
integer(pInt), intent(in) :: iFace
@ -1277,9 +1287,8 @@ function homogenization_RGC_grain1to3(&
)
use prec, only: pInt,p_vec
implicit none
!* Definition of variables
implicit none
integer(pInt), dimension (3) :: homogenization_RGC_grain1to3
integer(pInt), intent(in) :: grain1,homID
integer(pInt), dimension (3) :: nGDim
@ -1301,9 +1310,8 @@ function homogenization_RGC_grain3to1(&
)
use prec, only: pInt,p_vec
implicit none
!* Definition of variables
implicit none
integer(pInt), dimension (3), intent(in) :: grain3
integer(pInt) :: homogenization_RGC_grain3to1
integer(pInt), dimension (3) :: nGDim
@ -1324,9 +1332,8 @@ function homogenization_RGC_interface4to1(&
)
use prec, only: pInt,p_vec
implicit none
!* Definition of variables
implicit none
integer(pInt), dimension (4), intent(in) :: iFace4D
integer(pInt) :: homogenization_RGC_interface4to1
integer(pInt), dimension (3) :: nGDim,nIntFace
@ -1364,9 +1371,8 @@ function homogenization_RGC_interface1to4(&
)
use prec, only: pReal,pInt,p_vec
implicit none
!* Definition of variables
implicit none
integer(pInt), dimension (4) :: homogenization_RGC_interface1to4
integer(pInt), intent(in) :: iFace1D
integer(pInt), dimension (3) :: nGDim,nIntFace
@ -1442,9 +1448,8 @@ subroutine homogenization_RGC_grainDeformation(&
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance
implicit none
!* Definition of variables
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0
real(pReal), dimension (3,3), intent(in) :: avgF

86
code/homogenization_isostrain.f90
View File

@ -29,22 +29,27 @@
! Ngrains 6
! (output) Ngrains
MODULE homogenization_isostrain
module homogenization_isostrain
use prec, only: pInt
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
character (len=*), parameter :: &
homogenization_isostrain_label = 'isostrain'
character (len=*), parameter :: homogenization_isostrain_label = 'isostrain'
integer(pInt),dimension(:), allocatable :: &
homogenization_isostrain_sizeState, &
homogenization_isostrain_Ngrains, &
homogenization_isostrain_sizePostResults
integer(pInt), dimension(:), allocatable :: homogenization_isostrain_sizeState, &
homogenization_isostrain_Ngrains
integer(pInt), dimension(:), allocatable :: homogenization_isostrain_sizePostResults
integer(pInt), dimension(:,:), allocatable,target :: homogenization_isostrain_sizePostResult
character(len=64), dimension(:,:), allocatable,target :: homogenization_isostrain_output ! name of each post result output
integer(pInt), dimension(:,:), allocatable, target :: &
homogenization_isostrain_sizePostResult
character(len=64), dimension(:,:), allocatable, target :: &
homogenization_isostrain_output ! name of each post result output
CONTAINS
contains
!****************************************
!* - homogenization_isostrain_init
!* - homogenization_isostrain_stateInit
@ -58,9 +63,7 @@ CONTAINS
!**************************************
!* Module initialization *
!**************************************
subroutine homogenization_isostrain_init(&
myFile & ! file pointer to material configuration
)
subroutine homogenization_isostrain_init(myFile) ! file pointer to material configuration
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt
use math, only: math_Mandel3333to66, math_Voigt66to3333
@ -71,8 +74,8 @@ subroutine homogenization_isostrain_init(&
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
integer(pInt) section, i, j, output, mySize
integer :: maxNinstance, k !no pInt (stores a system dependen value from 'count'
character(len=64) tag
character(len=1024) line
character(len=64) :: tag
character(len=1024) :: line
!$OMP CRITICAL (write2out)
write(6,*)
@ -93,7 +96,6 @@ subroutine homogenization_isostrain_init(&
maxNinstance)) ; homogenization_isostrain_output = ''
rewind(myFile)
line = ''
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to <homogenization>
@ -144,9 +146,7 @@ subroutine homogenization_isostrain_init(&
enddo
enddo
return
endsubroutine
end subroutine homogenization_isostrain_init
!*********************************************************************
@ -154,18 +154,15 @@ endsubroutine
!*********************************************************************
function homogenization_isostrain_stateInit(myInstance)
use prec, only: pReal,pInt
implicit none
!* Definition of variables
implicit none
integer(pInt), intent(in) :: myInstance
real(pReal), dimension(homogenization_isostrain_sizeState(myInstance)) :: &
homogenization_isostrain_stateInit
homogenization_isostrain_stateInit = 0.0_pReal
return
endfunction
endfunction homogenization_isostrain_stateInit
!********************************************************************
@ -183,9 +180,8 @@ subroutine homogenization_isostrain_partitionDeformation(&
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains,homogenization_Ngrains
implicit none
!* Definition of variables
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0
real(pReal), dimension (3,3), intent(in) :: avgF
@ -197,9 +193,7 @@ subroutine homogenization_isostrain_partitionDeformation(&
forall (i = 1_pInt:homogenization_Ngrains(mesh_element(3,el))) &
F(1:3,1:3,i) = avgF
return
endsubroutine
end subroutine homogenization_isostrain_partitionDeformation
!********************************************************************
@ -230,9 +224,7 @@ function homogenization_isostrain_updateState(&
! homID = homogenization_typeInstance(mesh_element(3,el))
homogenization_isostrain_updateState = .true. ! homogenization at material point converged (done and happy)
return
endfunction
end function homogenization_isostrain_updateState
!********************************************************************
@ -251,9 +243,8 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(&
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains, homogenization_Ngrains
implicit none
!* Definition of variables
implicit none
real(pReal), dimension (3,3), intent(out) :: avgP
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P
@ -266,9 +257,7 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(&
avgP = sum(P,3)/real(Ngrains,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal)
return
endsubroutine
end subroutine homogenization_isostrain_averageStressAndItsTangent
!********************************************************************
@ -283,9 +272,8 @@ function homogenization_isostrain_averageTemperature(&
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains, homogenization_Ngrains
implicit none
!* Definition of variables
implicit none
real(pReal), dimension (homogenization_maxNgrains), intent(in) :: Temperature
integer(pInt), intent(in) :: ip,el
real(pReal) homogenization_isostrain_averageTemperature
@ -295,9 +283,7 @@ function homogenization_isostrain_averageTemperature(&
Ngrains = homogenization_Ngrains(mesh_element(3,el))
homogenization_isostrain_averageTemperature = sum(Temperature(1:Ngrains))/real(Ngrains,pReal)
return
endfunction
end function homogenization_isostrain_averageTemperature
!********************************************************************
@ -312,17 +298,15 @@ pure function homogenization_isostrain_postResults(&
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use material, only: homogenization_typeInstance,homogenization_Noutput
implicit none
!* Definition of variables
implicit none
type(p_vec), intent(in) :: state
integer(pInt), intent(in) :: ip,el
integer(pInt) homID,o,c
real(pReal), dimension(homogenization_isostrain_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: &
homogenization_isostrain_postResults
homID = homogenization_typeInstance(mesh_element(3,el))
integer(pInt) :: homID,o,c
real(pReal), dimension(homogenization_isostrain_sizePostResults&
(homogenization_typeInstance(mesh_element(3,el)))) :: homogenization_isostrain_postResults
c = 0_pInt
homID = homogenization_typeInstance(mesh_element(3,el))
homogenization_isostrain_postResults = 0.0_pReal
do o = 1_pInt,homogenization_Noutput(mesh_element(3,el))
@ -335,6 +319,6 @@ pure function homogenization_isostrain_postResults(&
return
endfunction
end function homogenization_isostrain_postResults
END MODULE
end module homogenization_isostrain

130
code/lattice.f90
View File

@ -27,56 +27,64 @@
!* - Schmid matrices calculation *
!************************************
MODULE lattice
module lattice
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
use prec, only: pReal,pInt
implicit none
!************************************
!* Lattice structures *
!************************************
integer(pInt) lattice_Nhexagonal, & ! # of hexagonal lattice structure (from tag CoverA_ratio)
lattice_Nstructure ! # of lattice structures (1: fcc,2: bcc,3+: hexagonal)
integer(pInt), parameter :: lattice_maxNslipFamily = 5_pInt ! max # of slip system families over lattice structures
integer(pInt), parameter :: lattice_maxNtwinFamily = 4_pInt ! max # of twin system families over lattice structures
integer(pInt), parameter :: lattice_maxNslip = 54_pInt ! max # of slip systems over lattice structures
integer(pInt), parameter :: lattice_maxNtwin = 24_pInt ! max # of twin systems over lattice structures
integer(pInt), parameter :: lattice_maxNinteraction = 30_pInt ! max # of interaction types (in hardening matrix part)
integer(pInt) :: &
lattice_Nhexagonal, & !> # of hexagonal lattice structure (from tag CoverA_ratio)
lattice_Nstructure !> # of lattice structures (1: fcc,2: bcc,3+: hexagonal)
integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, &
interactionSlipTwin, &
interactionTwinSlip, &
interactionTwinTwin
integer(pInt), parameter :: &
lattice_maxNslipFamily = 5_pInt, & !> max # of slip system families over lattice structures
lattice_maxNtwinFamily = 4_pInt, & !> max # of twin system families over lattice structures
lattice_maxNslip = 54_pInt, & !> max # of slip systems over lattice structures
lattice_maxNtwin = 24_pInt, & !> max # of twin systems over lattice structures
lattice_maxNinteraction = 30_pInt !> max # of interaction types (in hardening matrix part)
! Schmid matrices, normal, shear direction and d x n of slip systems
real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Sslip
real(pReal), allocatable, dimension(:,:,:) :: lattice_Sslip_v
real(pReal), allocatable, dimension(:,:,:) :: lattice_sn, &
lattice_sd, &
lattice_st
integer(pInt), pointer, dimension(:,:) :: &
interactionSlipSlip, &
interactionSlipTwin, &
interactionTwinSlip, &
interactionTwinTwin
real(pReal), allocatable, dimension(:,:,:,:) :: &
lattice_Sslip ! Schmid matrices, normal, shear direction and d x n of slip systems
real(pReal), allocatable, dimension(:,:,:) :: &
lattice_Sslip_v, &
lattice_sn, &
lattice_sd, &
lattice_st
! rotation and Schmid matrices, normal, shear direction and d x n of twin systems
real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Qtwin
real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Stwin
real(pReal), allocatable, dimension(:,:,:) :: lattice_Stwin_v
real(pReal), allocatable, dimension(:,:,:) :: lattice_tn, &
lattice_td, &
lattice_tt
real(pReal), allocatable, dimension(:,:,:,:) :: &
lattice_Qtwin, &
lattice_Stwin
! characteristic twin shear
real(pReal), allocatable, dimension(:,:) :: lattice_shearTwin
real(pReal), allocatable, dimension(:,:,:) :: &
lattice_Stwin_v, &
lattice_tn, &
lattice_td, &
lattice_tt
! number of slip and twin systems in each family
integer(pInt), allocatable, dimension(:,:) :: lattice_NslipSystem, &
lattice_NtwinSystem
real(pReal), allocatable, dimension(:,:) :: &
lattice_shearTwin !> characteristic twin shear
! interaction type of slip and twin systems among each other
integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip, &
lattice_interactionSlipTwin, &
lattice_interactionTwinSlip, &
lattice_interactionTwinTwin
integer(pInt), allocatable, dimension(:,:) :: &
lattice_NslipSystem, & !> number of slip systems in each family
lattice_NtwinSystem !> number of twin systems in each family
integer(pInt), allocatable, dimension(:,:,:) :: &
lattice_interactionSlipSlip, & !> interaction type between slip/slip
lattice_interactionSlipTwin, & !> interaction type between slip/twin
lattice_interactionTwinSlip, & !> interaction type between twin/slip
lattice_interactionTwinTwin !> interaction type between twin/twin
!============================== fcc (1) =================================
@ -698,16 +706,15 @@ CONTAINS
!* - lattice_initializeStructure
!****************************************
pure function lattice_symmetryType(structID)
integer(pInt) pure function lattice_symmetryType(structID)
!**************************************
!* maps structure to symmetry type *
!* fcc(1) and bcc(2) are cubic(1) *
!* hex(3+) is hexagonal(2) *
!**************************************
implicit none
implicit none
integer(pInt), intent(in) :: structID
integer(pInt) lattice_symmetryType
select case(structID)
case (1_pInt,2_pInt)
@ -720,21 +727,29 @@ pure function lattice_symmetryType(structID)
return
end function
end function lattice_symmetryType
subroutine lattice_init()
subroutine lattice_init
!**************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use IO, only: IO_open_file,IO_open_jobFile_stat,IO_countSections,IO_countTagInPart,IO_error
use material, only: material_configfile,material_localFileExt,material_partPhase
use debug, only: debug_verbosity
implicit none
use IO, only: IO_open_file,&
IO_open_jobFile_stat, &
IO_countSections, &
IO_countTagInPart, &
IO_error
use material, only: material_configfile, &
material_localFileExt, &
material_partPhase
use debug, only: debug_what, &
debug_lattice, &
debug_levelBasic
implicit none
integer(pInt), parameter :: fileunit = 200_pInt
integer(pInt) Nsections
integer(pInt) :: Nsections
!$OMP CRITICAL (write2out)
write(6,*)
@ -751,7 +766,7 @@ subroutine lattice_init()
! lattice_Nstructure = Nsections + 2_pInt ! most conservative assumption
close(fileunit)
if (debug_verbosity > 0_pInt) then
if (iand(debug_what(debug_lattice),debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a16,1x,i5)') '# phases:',Nsections
write(6,'(a16,1x,i5)') '# structures:',lattice_Nstructure
@ -782,19 +797,25 @@ subroutine lattice_init()
allocate(lattice_interactionTwinSlip(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt ! other:me
allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt ! other:me
end subroutine
end subroutine lattice_init
function lattice_initializeStructure(struct,CoverA)
integer(pInt) function lattice_initializeStructure(struct,CoverA)
!**************************************
!* Calculation of Schmid *
!* matrices, etc. *
!**************************************
use prec, only: pReal,pInt
use math
use math, only: math_vectorproduct, &
math_tensorproduct, &
math_mul3x3, &
math_symmetric33, &
math_Mandel33to6, &
math_axisAngleToR, &
INRAD
use IO, only: IO_error
implicit none
implicit none
character(len=*) struct
real(pReal) CoverA
real(pReal), dimension(3,lattice_maxNslip) :: sd = 0.0_pReal, &
@ -811,7 +832,6 @@ function lattice_initializeStructure(struct,CoverA)
integer(pInt) :: i,myNslip,myNtwin,myStructure = 0_pInt
logical :: processMe
integer(pInt) lattice_initializeStructure
processMe = .false.
select case(struct(1:3)) ! check first three chars of structure name
@ -949,7 +969,7 @@ function lattice_initializeStructure(struct,CoverA)
lattice_initializeStructure = myStructure ! report my structure index back
end function
end function lattice_initializeStructure
END MODULE
end module lattice

355
code/makefile
View File

@ -1,355 +0,0 @@
########################################################################################
# Makefile to compile the Material subroutine for BVP solution using spectral method
########################################################################################
# Be sure to remove all files compiled with different options by using "make clean"
#
# Uses OpenMP to parallelize the material subroutines (set number of threads with "export DAMASK_NUM_THREADS=n" to n)
#
# Install fftw3 (v3.3 is tested):
# + run
# ./configure --enable-threads --enable-sse2 --enable-shared [-enable-float]
# make
# make install
# + specify in the "pathinfo:FFTW" where FFTW was installed.
# We essentially look for two library files "lib/libfftw3_threads" and "lib/libfftw3", so you can copy those, for instance,
# into DAMASK_ROOT/lib/fftw/lib/ and specify "./fftw/" as pathinfo:FFTW
# Use --enable-float in above configure for single precision...
# Uses linux threads to parallelize fftw3
#
# Instead of the AMD Core Math Library a standard "liblapack.a/dylib/etc." can be used by leaving pathinfo:ACML and pathinfo:IKML blank
########################################################################################
# OPTIONS = standard (alternative): meaning
#-------------------------------------------------------------
# F90 = ifort (gfortran): compiler, choose Intel or GNU
# COMPILERNAME = overwrite name of Compiler, e.g. using mpich-g90 instead of ifort
# PORTABLE = TRUE (FALSE): decision, if executable is optimized for the machine on which it was built.
# OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE,ULTRA): Optimization mode: O2, O0, O3 + further options for most files, O3 + further options for all files
# OPENMP = TRUE (FALSE): OpenMP multiprocessor support
# FFTWROOT = pathinfo:FFTW (will be adjusted by setup_code.py - required in pathinfo)
# IKMLROOT = pathinfo:IKML (will be adjusted by setup_code.py if present in pathinfo)
# ACMLROOT = pathinfo:ACML (will be adjusted by setup_code.py if present in pathinfo)
# LAPACKROOT = pathinfo:LAPACK (will be adjusted by setup_code.py if present in pathinfo)
# PREFIX = arbitrary prefix
# SUFFIX = arbitrary suffix
# STANDARD_CHECK = checking for Fortran 2008, compiler dependend
########################################################################################
#auto values will be set by setup_code.py
FFTWROOT :=/$(DAMASK_ROOT)/lib/fftw
IKMLROOT :=
ACMLROOT :=/opt/acml4.4.0
#LAPACKROOT := /usr
F90 ?= ifort
COMPILERNAME ?= $(F90)
OPENMP ?= ON
OPTIMIZATION ?= DEFENSIVE
ifeq "$(F90)" "ifort"
ifeq "$(OPTIMIZATION)" "OFF"
ARCHIVE_COMMAND :=ar
else
ARCHIVE_COMMAND :=xiar
endif
else
ARCHIVE_COMMAND :=ar
endif
ifeq "$(OPTIMIZATION)" "OFF"
OPTI := OFF
MAXOPTI := OFF
endif
ifeq "$(OPTIMIZATION)" "DEFENSIVE"
OPTI := DEFENSIVE
MAXOPTI := DEFENSIVE
endif
ifeq "$(OPTIMIZATION)" "AGGRESSIVE"
OPTI := AGGRESSIVE
MAXOPTI := DEFENSIVE
endif
ifeq "$(OPTIMIZATION)" "ULTRA"
OPTI := AGGRESSIVE
MAXOPTI := AGGRESSIVE
endif
ifndef OPTI
OPTI := DEFENSIVE
MAXOPTI := DEFENSIVE
endif
ifeq "$(PORTABLE)" "FALSE"
PORTABLE_SWITCH =-msse3
endif
# settings for multicore support
ifeq "$(OPENMP)" "ON"
OPENMP_FLAG_ifort =-openmp -openmp-report0 -parallel
OPENMP_FLAG_gfortran =-fopenmp
ACML_ARCH =_mp
LIBRARIES +=-lfftw3_threads -lpthread
endif
LIBRARIES +=-lfftw3
LIB_DIRS +=-L$(FFTWROOT)/lib
ifdef IKMLROOT
LIBRARIES +=-mkl
else
ifdef ACMLROOT
LIB_DIRS +=-L$(ACMLROOT)/$(F90)64$(ACML_ARCH)/lib
LIBRARIES +=-lacml$(ACML_ARCH)
else
ifdef LAPACKROOT
LIB_DIRS +=-L$(LAPACKROOT)/lib64 -L$(LAPACKROOT)/lib
LIBRARIES +=-llapack
endif
endif
endif
ifdef STANDARD_CHECK
STANDARD_CHECK_ifort =$(STANDARD_CHECK)
STANDARD_CHECK_gfortran =$(STANDARD_CHECK)
endif
STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics
STANDARD_CHECK_gfortran ?=-std=f2008
OPTIMIZATION_OFF_ifort :=-O0
OPTIMIZATION_OFF_gfortran :=-O0
OPTIMIZATION_DEFENSIVE_ifort :=-O2
OPTIMIZATION_DEFENSIVE_gfortran :=-O2
OPTIMIZATION_AGGRESSIVE_ifort :=-O3 $(PORTABLE_SWITCH) -ipo -static -no-prec-div -fp-model fast=2
OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 $(PORTABLE_SWITCH) -ffast-math -funroll-loops -ftree-vectorize
COMPILE_OPTIONS_ifort :=-fpp\
-implicitnone\
-diag-enable sc3\
-diag-disable 5268\
-warn declarations\
-warn general\
-warn usage\
-warn interfaces\
-warn ignore_loc\
-warn alignments\
-warn unused\
-warn errors\
-warn stderrors
#-fpp: preprocessor
#-fimplicit-none: assume "implicit-none" even if not present in source
#-diag-disable: disables warnings, where
# warning ID 5268: the text exceeds right hand column allowed on the line (we have only comments there)
#-warn: enables warnings, where
# declarations: any undeclared names
# general: warning messages and informational messages are issued by the compiler
# usage: questionable programming practices
# interfaces: checks the interfaces of all SUBROUTINEs called and FUNCTIONs invoked in your compilation against an external set of interface blocks
# ignore_loc: %LOC is stripped from an actual argument
# alignments: data that is not naturally aligned
# unused: declared variables that are never used
# errors: warnings are changed to errors
# stderrors: warnings about Fortran standard violations are changed to errors
#
###################################################################################################
#MORE OPTIONS FOR DEBUGGING DURING COMPILING
#-warn: enables warnings, where
# truncated_source: Determines whether warnings occur when source exceeds the maximum column width in fixed-format files. (too many warnings because we have comments beyond character 132)
# uncalled: Determines whether warnings occur when a statement function is never called
# all:
#
#OPTIONS FOR DEGUBBING DURING RUNTIME
# information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/
#-g: Generate symbolic debugging information in the object file
#-traceback: Generate extra information in the object file to provide source file traceback information when a severe error occurs at run time.
#-gen-interfaces: Generate an interface block for each routine. http://software.intel.com/en-us/blogs/2012/01/05/doctor-fortran-gets-explicit-again/
#-fp-stack-check: Generate extra code after every function call to ensure that the floating-point (FP) stack is in the expected state.
#-check: checks at runtime, where
# bounds: check if an array index is too small (<1) or too large!
# arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays
# format: Checking for the data type of an item being formatted for output.
# output_conversion: Checking for the fit of data items within a designated format descriptor field.
# pointers: Checking for certain disassociated or uninitialized pointers or unallocated allocatable objects.
# uninit: Checking for uninitialized variables.
#-heap-arrays: should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
#
#OPTIONS FOR TYPE DEBUGGING
#-real-size 32: set precision to one of those 32/64/128 (= 4/8/16 bytes) for standard real (=8 for pReal)
#-integer-size 16: set precision to one of those 16/32/64 (= 2/4/8 bytes) for standard integer (=4 for pInt)
###################################################################################################
COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\
-ffree-line-length-132\
-fno-range-check\
-fimplicit-none\
-fall-intrinsics\
-pedantic\
-Warray-bounds\
-Wampersand\
-Wno-tabs\
-Wcharacter-truncation\
-Wintrinsic-shadow\
-Waliasing\
-Wconversion\
-Wsurprising\
-Wunderflow\
-Wswitch\
-Wstrict-overflow\
-Wattributes\
-Wunsafe-loop-optimizations\
-Wunused\
-Wextra
#-xf95-cpp-input: preprocessor
#-ffree-line-length-132: restrict line length to the standard 132 characters
#-fno-range-check: disables checking if result can be represented by variable. Needs to be set to enable DAMASK_NaN
#-fimplicit-none: assume "implicit-none" even if not present in source
#-fall-intrinsics:
#-pedantic: more strict on standard, enables some of the warnings below
#-Warray-bounds: checks if array reference is out of bounds at compile time. use -fcheck-bounds to also check during runtime
#-Wampersand: checks if a character expression is continued proberly by an ampersand at the end of the line and at the beginning of the new line
#-Wno-tabs: do not allow tabs in source
#-Wcharacter-truncation: warn if character expressions (strings) are truncated
#-Wintrinsic-shadow: warn if a user-defined procedure or module procedure has the same name as an intrinsic
#-Waliasing: warn about possible aliasing of dummy arguments. Specifically, it warns if the same actual argument is associated with a dummy argument with "INTENT(IN)" and a dummy argument with "INTENT(OUT)" in a call with an explicit interface.
#-Wconversion: warn about implicit conversions between different type
#-Wsurprising: warn when "suspicious" code constructs are encountered. While technically legal these usually indicate that an error has been made.
#-Wunderflow: produce a warning when numerical constant expressions are encountered, which yield an UNDERFLOW during compilation
#-Wswitch: warn whenever a "switch" statement has an index of enumerated type and lacks a "case" for one or more of the named codes of that enumeration. (The presence of a "default" label prevents this warning.) "case" labels outside the enumeration range also provokewarnings when this option is used (even if there is a "default" label)
#-Wstrict-overflow:
#-Wattributes: warn about inappropriate attribute usage
#-Wunsafe-loop-optimizations: warn if the loop cannot be optimized due to nontrivial assumptions.
#-Wunused:
# -value:
# -parameter: find usused variables with "parameter" attribute
#-Wextra:
###################################################################################################
#OPTIONS FOR GFORTRAN 4.6
#-Wsuggest-attribute=const:
#-Wsuggest-attribute=noreturn:
#-Wsuggest-attribute=pure:
#-Wreal-q-constant: Warn about real-literal-constants with 'q' exponent-letter
#MORE OPTIONS FOR DEBUGGING DURING COMPILING
#-Wline-truncation: too many warnings because we have comments beyond character 132
#-Wintrinsic-std: warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers:
#-Warray-temporarieswarnings:
# because we have many temporary arrays (performance issue?):
#-Wimplicit-interface
#-pedantic-errors
#-fmodule-private
#
#OPTIONS FOR DEGUBBING DURING RUNTIME
#-fcheck-bounds: check if an array index is too small (<1) or too large!
#
#OPTIONS FOR TYPE DEBUGGING
#-fdefault-real-8: set precision to 8 bytes for standard real (=8 for pReal). Will set size of double to 16 bytes as long as -fdefault-double-8 is not set
#-fdefault-integer-8: set precision to 8 bytes for standard integer (=4 for pInt)
##################################################################################################
COMPILE =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) -c
COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) -c
###################################################################################################
DAMASK_spectral.exe: DAMASK_spectral.o CPFEM.a
$(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a \
constitutive.a advanced.a basics.a $(LIB_DIRS) $(LIBRARIES)
DAMASK_spectral.o: DAMASK_spectral.f90 CPFEM.o
$(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) DAMASK_spectral.f90 $(SUFFIX)
CPFEM.a: CPFEM.o
$(ARCHIVE_COMMAND) rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o
CPFEM.o: CPFEM.f90 homogenization.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) CPFEM.f90 $(SUFFIX)
homogenization.o: homogenization.f90 homogenization_isostrain.o homogenization_RGC.o crystallite.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) homogenization.f90 $(SUFFIX)
homogenization_RGC.o: homogenization_RGC.f90 constitutive.a
$(PREFIX) $(COMPILERNAME) $(COMPILE) homogenization_RGC.f90 $(SUFFIX)
homogenization_isostrain.o: homogenization_isostrain.f90 basics.a advanced.a
$(PREFIX) $(COMPILERNAME) $(COMPILE) homogenization_isostrain.f90 $(SUFFIX)
crystallite.o: crystallite.f90 constitutive.a
$(PREFIX) $(COMPILERNAME) $(COMPILE) crystallite.f90 $(SUFFIX)
constitutive.a: constitutive.o
$(ARCHIVE_COMMAND) rc constitutive.a constitutive.o constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o basics.a advanced.a
constitutive.o: constitutive.f90 constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive.f90 $(SUFFIX)
constitutive_titanmod.o: constitutive_titanmod.f90 basics.a advanced.a
$(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_titanmod.f90 $(SUFFIX)
constitutive_nonlocal.o: constitutive_nonlocal.f90 basics.a advanced.a
$(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_nonlocal.f90 $(SUFFIX)
constitutive_dislotwin.o: constitutive_dislotwin.f90 basics.a advanced.a
$(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_dislotwin.f90 $(SUFFIX)
constitutive_j2.o: constitutive_j2.f90 basics.a advanced.a
$(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_j2.f90 $(SUFFIX)
constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 basics.a advanced.a
$(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_phenopowerlaw.f90 $(SUFFIX)
advanced.a: lattice.o
$(ARCHIVE_COMMAND) rc advanced.a FEsolving.o mesh.o material.o lattice.o
lattice.o: lattice.f90 material.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) lattice.f90 $(SUFFIX)
material.o: material.f90 mesh.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) material.f90 $(SUFFIX)
mesh.o: mesh.f90 FEsolving.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) mesh.f90 $(SUFFIX)
FEsolving.o: FEsolving.f90 basics.a
$(PREFIX) $(COMPILERNAME) $(COMPILE) FEsolving.f90 $(SUFFIX)
basics.a: math.o
$(ARCHIVE_COMMAND) rc basics.a math.o debug.o numerics.o IO.o DAMASK_spectral_interface.o prec.o
math.o: math.f90 debug.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) math.f90 $(SUFFIX)
debug.o: debug.f90 numerics.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) debug.f90 $(SUFFIX)
numerics.o: numerics.f90 IO.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) numerics.f90 $(SUFFIX)
IO.o: IO.f90 DAMASK_spectral_interface.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) IO.f90 $(SUFFIX)
DAMASK_spectral_interface.o: DAMASK_spectral_interface.f90 prec.o
$(PREFIX) $(COMPILERNAME) $(COMPILE) DAMASK_spectral_interface.f90 $(SUFFIX)
prec.o: prec.f90
$(PREFIX) $(COMPILERNAME) $(COMPILE) prec.f90 $(SUFFIX)
tidy:
rm -rf *.o
rm -rf *.mod
rm -rf *.a
clean:
rm -rf *.o
rm -rf *.mod
rm -rf *.a
rm -rf *.exe

341
code/material.f90
View File

@ -16,106 +16,129 @@
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!##############################################################
!--------------------------------------------------------------------------------------------------
!* $Id$
!************************************
!* Module: MATERIAL *
!************************************
!* contains: *
!* - parsing of material.config *
!************************************
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!! Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Parses material.config
!--------------------------------------------------------------------------------------------------
module material
MODULE material
use prec, only: pReal, &
pInt
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
implicit none
private
character(len=64), parameter, public :: &
material_configFile = 'material.config', &
material_localFileExt = 'materialConfig'
character(len=64), parameter, public :: material_configFile = 'material.config'
character(len=64), parameter, public :: material_localFileExt = 'materialConfig'
character(len=32), parameter, public :: material_partHomogenization = 'homogenization'
character(len=32), parameter, private :: material_partMicrostructure = 'microstructure'
character(len=32), parameter, public :: material_partCrystallite = 'crystallite'
character(len=32), parameter, public :: material_partPhase = 'phase'
character(len=32), parameter, private :: material_partTexture = 'texture'
character(len=32), parameter, public :: &
material_partHomogenization = 'homogenization', &
material_partCrystallite = 'crystallite', &
material_partPhase = 'phase'
character(len=64), dimension(:), allocatable, public :: &
phase_constitution, & !> constitution of each phase
phase_name, & !> name of each phase
homogenization_name, & !> name of each homogenization
homogenization_type, & !> type of each homogenization
crystallite_name !> name of each crystallite setting
integer(pInt), public :: &
homogenization_maxNgrains, & !> max number of grains in any USED homogenization
material_Nphase, & !> number of phases
material_Nhomogenization, & !> number of homogenizations
material_Nmicrostructure, & !> number of microstructures
material_Ncrystallite !> number of crystallite settings
integer(pInt), dimension(:), allocatable, public :: &
homogenization_Ngrains, & !> number of grains in each homogenization
homogenization_Noutput, & !> number of '(output)' items per homogenization
phase_Noutput, & !> number of '(output)' items per phase
phase_constitutionInstance, & !> instance of particular constitution of each phase
crystallite_Noutput, & !> number of '(output)' items per crystallite setting
homogenization_typeInstance, & !> instance of particular type of each homogenization
microstructure_crystallite !> crystallite setting ID of each microstructure
integer(pInt), dimension(:,:,:), allocatable, public :: &
material_phase, & !> phase (index) of each grain,IP,element
material_texture !> texture (index) of each grain,IP,element
real(pReal), dimension(:,:,:,:), allocatable, public :: &
material_EulerAngles !> initial orientation of each grain,IP,element
logical, dimension(:), allocatable, public :: &
microstructure_active, &
microstructure_elemhomo, & !> flag to indicate homogeneous microstructure distribution over element's IPs
phase_localConstitution !> flags phases with local constitutive law
!*************************************
!* Definition of material properties *
!*************************************
!* Number of materials
integer(pInt) &
material_Nhomogenization, & ! number of homogenizations
material_Nmicrostructure, & ! number of microstructures
material_Ncrystallite, & ! number of crystallite settings
material_Nphase, & ! number of phases
material_Ntexture, & ! number of textures
microstructure_maxNconstituents,&! max number of constituents in any phase
homogenization_maxNgrains, & ! max number of grains in any USED homogenization
texture_maxNgauss, & ! max number of Gauss components in any texture
texture_maxNfiber ! max number of Fiber components in any texture
character(len=64), dimension(:), allocatable :: &
homogenization_name, & ! name of each homogenization
homogenization_type, & ! type of each homogenization
microstructure_name, & ! name of each microstructure
crystallite_name, & ! name of each crystallite setting
phase_name, & ! name of each phase
phase_constitution, & ! constitution of each phase
texture_name ! name of each texture
character(len=256),dimension(:), allocatable :: &
texture_ODFfile ! name of each ODF file
integer(pInt), dimension(:), allocatable :: &
homogenization_Ngrains, & ! number of grains in each homogenization
homogenization_typeInstance, & ! instance of particular type of each homogenization
homogenization_Noutput, & ! number of '(output)' items per homogenization
microstructure_Nconstituents, & ! number of constituents in each microstructure
crystallite_Noutput, & ! number of '(output)' items per crystallite setting
phase_constitutionInstance, & ! instance of particular constitution of each phase
phase_Noutput, & ! number of '(output)' items per phase
texture_symmetry, & ! number of symmetric orientations per texture
texture_Ngauss, & ! number of Gauss components per texture
texture_Nfiber ! number of Fiber components per texture
logical, dimension(:), allocatable :: &
homogenization_active, & !
microstructure_active, & !
microstructure_elemhomo, & ! flag to indicate homogeneous microstructure distribution over element's IPs
phase_localConstitution ! flags phases with local constitutive law
integer(pInt), dimension(:), allocatable :: &
microstructure_crystallite ! crystallite setting ID of each microstructure
integer(pInt), dimension(:,:), allocatable :: &
microstructure_phase, & ! phase IDs of each microstructure
microstructure_texture ! texture IDs of each microstructure
real(pReal), dimension(:,:), allocatable :: &
microstructure_fraction ! vol fraction of each constituent in microstructure
real(pReal), dimension(:,:,:), allocatable :: &
material_volume ! volume of each grain,IP,element
integer(pInt), dimension(:,:,:), allocatable :: &
material_phase, & ! phase (index) of each grain,IP,element
material_texture ! texture (index) of each grain,IP,element
real(pReal), dimension(:,:,:,:), allocatable :: &
material_EulerAngles ! initial orientation of each grain,IP,element
real(pReal), dimension(:,:,:), allocatable :: &
texture_Gauss, & ! data of each Gauss component
texture_Fiber ! data of each Fiber component
character(len=32), parameter, private :: &
material_partMicrostructure = 'microstructure', &
material_partTexture = 'texture'
CONTAINS
character(len=64), dimension(:), allocatable, private :: &
microstructure_name, & !> name of each microstructure
texture_name !> name of each texture
character(len=256), dimension(:), allocatable, private :: &
texture_ODFfile !> name of each ODF file
integer(pInt), private :: &
material_Ntexture, & !> number of textures
microstructure_maxNconstituents, & !> max number of constituents in any phase
texture_maxNgauss, & !> max number of Gauss components in any texture
texture_maxNfiber !> max number of Fiber components in any texture
integer(pInt), dimension(:), allocatable, private :: &
microstructure_Nconstituents, & !> number of constituents in each microstructure
texture_symmetry, & !> number of symmetric orientations per texture
texture_Ngauss, & !> number of Gauss components per texture
texture_Nfiber !> number of Fiber components per texture
integer(pInt), dimension(:,:), allocatable, private :: &
microstructure_phase, & !> phase IDs of each microstructure
microstructure_texture !> texture IDs of each microstructure
real(pReal), dimension(:,:), allocatable, private :: &
microstructure_fraction !> vol fraction of each constituent in microstructure
real(pReal), dimension(:,:,:), allocatable :: &
material_volume, & !> volume of each grain,IP,element
texture_Gauss, & !> data of each Gauss component
texture_Fiber !> data of each Fiber component
logical, dimension(:), allocatable, private :: &
homogenization_active
public :: material_init
contains
!*********************************************************************
subroutine material_init()
subroutine material_init
!*********************************************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pReal,pInt
use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat
use debug, only: debug_verbosity
implicit none
use IO, only: IO_error, &
IO_open_file, &
IO_open_jobFile_stat
use debug, only: debug_what, &
debug_material, &
debug_levelBasic, &
debug_levelExtensive
implicit none
!* Definition of variables
integer(pInt), parameter :: fileunit = 200_pInt
integer(pInt) i,j
integer(pInt) :: i,j, myDebug
myDebug = debug_what(debug_material)
!$OMP CRITICAL (write2out)
write(6,*)
@ -128,31 +151,31 @@ subroutine material_init()
call IO_open_file(fileunit,material_configFile) ! ...open material.config file
endif
call material_parseHomogenization(fileunit,material_partHomogenization)
if (debug_verbosity > 0_pInt) then
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write (6,*) 'Homogenization parsed'
!$OMP END CRITICAL (write2out)
endif
call material_parseMicrostructure(fileunit,material_partMicrostructure)
if (debug_verbosity > 0_pInt) then
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write (6,*) 'Microstructure parsed'
!$OMP END CRITICAL (write2out)
endif
call material_parseCrystallite(fileunit,material_partCrystallite)
if (debug_verbosity > 0_pInt) then
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write (6,*) 'Crystallite parsed'
!$OMP END CRITICAL (write2out)
endif
call material_parseTexture(fileunit,material_partTexture)
if (debug_verbosity > 0_pInt) then
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write (6,*) 'Texture parsed'
!$OMP END CRITICAL (write2out)
endif
call material_parsePhase(fileunit,material_partPhase)
if (debug_verbosity > 0_pInt) then
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write (6,*) 'Phase parsed'
!$OMP END CRITICAL (write2out)
@ -167,7 +190,7 @@ subroutine material_init()
if (minval(microstructure_texture(1:microstructure_Nconstituents(i),i)) < 1_pInt .or. &
maxval(microstructure_texture(1:microstructure_Nconstituents(i),i)) > material_Ntexture) call IO_error(152_pInt,i)
if (abs(sum(microstructure_fraction(:,i)) - 1.0_pReal) >= 1.0e-10_pReal) then
if (debug_verbosity > 0_pInt) then
if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,*)'sum of microstructure fraction = ',sum(microstructure_fraction(:,i))
!$OMP END CRITICAL (write2out)
@ -175,7 +198,7 @@ subroutine material_init()
call IO_error(153_pInt,i)
endif
enddo
if (debug_verbosity > 0_pInt) then
if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,*) 'MATERIAL configuration'
@ -188,14 +211,14 @@ subroutine material_init()
write (6,'(a32,1x,a11,1x,a12,1x,a13)') 'microstructure ','crystallite','constituents','homogeneous'
do i = 1_pInt,material_Nmicrostructure
write (6,'(a32,4x,i4,8x,i4,8x,l1)') microstructure_name(i), &
microstructure_crystallite(i), &
microstructure_Nconstituents(i), &
microstructure_elemhomo(i)
microstructure_crystallite(i), &
microstructure_Nconstituents(i), &
microstructure_elemhomo(i)
if (microstructure_Nconstituents(i) > 0_pInt) then
do j = 1_pInt,microstructure_Nconstituents(i)
write (6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(j,i)),&
texture_name(microstructure_texture(j,i)),&
microstructure_fraction(j,i)
texture_name(microstructure_texture(j,i)),&
microstructure_fraction(j,i)
enddo
write (6,*)
endif
@ -203,27 +226,28 @@ subroutine material_init()
!$OMP END CRITICAL (write2out)
endif
call material_populateGrains()
call material_populateGrains
endsubroutine
end subroutine material_init
!*********************************************************************
subroutine material_parseHomogenization(myFile,myPart)
!*********************************************************************
use prec, only: pInt
use IO
use mesh, only: mesh_element
implicit none
implicit none
character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 2_pInt
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 2_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section, s
character(len=64) tag
character(len=1024) line
character(len=64) :: tag
character(len=1024) ::line
Nsections = IO_countSections(myFile,myPart)
material_Nhomogenization = Nsections
@ -273,25 +297,26 @@ subroutine material_parseHomogenization(myFile,myPart)
100 homogenization_maxNgrains = maxval(homogenization_Ngrains,homogenization_active)
endsubroutine
end subroutine material_parseHomogenization
!*********************************************************************
subroutine material_parseMicrostructure(myFile,myPart)
!*********************************************************************
use prec, only: pInt
use IO
use mesh, only: mesh_element, mesh_NcpElems
implicit none
implicit none
character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: myFile
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 7_pInt
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
integer(pInt) Nsections, section, constituent, e, i
character(len=64) tag
character(len=1024) line
integer(pInt) :: Nsections, section, constituent, e, i
character(len=64) :: tag
character(len=1024) :: line
Nsections = IO_countSections(myFile,myPart)
material_Nmicrostructure = Nsections
@ -353,21 +378,27 @@ subroutine material_parseMicrostructure(myFile,myPart)
endif
enddo
100 endsubroutine
100 end subroutine material_parseMicrostructure
!*********************************************************************
subroutine material_parseCrystallite(myFile,myPart)
!*********************************************************************
use prec, only: pInt
use IO
implicit none
use IO, only: IO_countSections, &
IO_error, &
IO_countTagInPart, &
IO_getTag, &
IO_lc, &
IO_isBlank
implicit none
character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: myFile
integer(pInt) Nsections, section
character(len=1024) line
integer(pInt), intent(in) :: myFile
integer(pInt) :: Nsections, &
section
character(len=1024) :: line
Nsections = IO_countSections(myFile,myPart)
material_Ncrystallite = Nsections
@ -396,24 +427,25 @@ subroutine material_parseCrystallite(myFile,myPart)
endif
enddo
100 endsubroutine
100 end subroutine material_parseCrystallite
!*********************************************************************
subroutine material_parsePhase(myFile,myPart)
!*********************************************************************
use prec, only: pInt
use IO
implicit none
implicit none
character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: myFile
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 2_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section, s
character(len=64) tag
character(len=1024) line
character(len=64) :: tag
character(len=1024) :: line
Nsections = IO_countSections(myFile,myPart)
material_Nphase = Nsections
@ -458,25 +490,26 @@ subroutine material_parsePhase(myFile,myPart)
endif
enddo
100 endsubroutine
100 end subroutine material_parsePhase
!*********************************************************************
subroutine material_parseTexture(myFile,myPart)
!*********************************************************************
use prec, only: pInt, pReal
use IO
use math, only: inRad, math_sampleRandomOri
implicit none
implicit none
character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 13_pInt
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 13_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section, gauss, fiber, i
character(len=64) tag
character(len=1024) line
integer(pInt) :: Nsections, section, gauss, fiber, i
character(len=64) :: tag
character(len=1024) :: line
Nsections = IO_countSections(myFile,myPart)
@ -589,35 +622,46 @@ subroutine material_parseTexture(myFile,myPart)
endif
enddo
100 endsubroutine
100 end subroutine material_parseTexture
!*********************************************************************
subroutine material_populateGrains()
subroutine material_populateGrains
!*********************************************************************
use prec, only: pInt, pReal
use math, only: math_sampleRandomOri, math_sampleGaussOri, math_sampleFiberOri, math_symmetricEulers
use mesh, only: mesh_element, mesh_maxNips, mesh_NcpElems, mesh_ipVolume, FE_Nips
use IO, only: IO_error, IO_hybridIA
use math, only: math_sampleRandomOri, &
math_sampleGaussOri, &
math_sampleFiberOri, &
math_symmetricEulers
use mesh, only: mesh_element, &
mesh_maxNips, &
mesh_NcpElems, &
mesh_ipVolume, &
FE_Nips
use IO, only: IO_error, &
IO_hybridIA
use FEsolving, only: FEsolving_execIP
use debug, only: debug_verbosity
implicit none
use debug, only: debug_what, &
debug_material, &
debug_levelBasic
implicit none
integer(pInt), dimension (:,:), allocatable :: Ngrains
integer(pInt), dimension (microstructure_maxNconstituents) :: NgrainsOfConstituent
integer(pInt), dimension (microstructure_maxNconstituents) &
:: NgrainsOfConstituent
real(pReal), dimension (:), allocatable :: volumeOfGrain
real(pReal), dimension (:,:), allocatable :: orientationOfGrain
real(pReal), dimension (3) :: orientation
real(pReal), dimension (3,3) :: symOrientation
real(pReal), dimension (3) :: orientation
real(pReal), dimension (3,3) :: symOrientation
integer(pInt), dimension (:), allocatable :: phaseOfGrain, textureOfGrain
integer(pInt) t,e,i,g,j,m,homog,micro,sgn,hme
integer(pInt) phaseID,textureID,dGrains,myNgrains,myNorientations, &
integer(pInt) :: t,e,i,g,j,m,homog,micro,sgn,hme, myDebug
integer(pInt) :: phaseID,textureID,dGrains,myNgrains,myNorientations, &
grain,constituentGrain,symExtension
real(pReal) extreme,rnd
real(pReal) :: extreme,rnd
integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array
integer(pInt), dimension (:,:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array
myDebug = debug_what(debug_material)
allocate(material_volume(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; material_volume = 0.0_pReal
allocate(material_phase(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; material_phase = 0_pInt
@ -663,7 +707,7 @@ subroutine material_populateGrains()
allocate(textureOfGrain(maxval(Ngrains))) ! reserve memory for maximum case
allocate(orientationOfGrain(3,maxval(Ngrains))) ! reserve memory for maximum case
if (debug_verbosity > 0_pInt) then
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,*) 'MATERIAL grain population'
@ -676,7 +720,7 @@ subroutine material_populateGrains()
do micro = 1_pInt,material_Nmicrostructure ! all pairs of homog and micro
if (Ngrains(homog,micro) > 0_pInt) then ! an active pair of homog and micro
myNgrains = Ngrains(homog,micro) ! assign short name for total number of grains to populate
if (debug_verbosity > 0_pInt) then
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a32,1x,a32,1x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains
@ -837,7 +881,6 @@ subroutine material_populateGrains()
deallocate(Nelems)
deallocate(elemsOfHomogMicro)
endsubroutine
end subroutine material_populateGrains
END MODULE
end module material

521
code/math.f90
View File

File diff suppressed because it is too large Load Diff

1197
code/mesh.f90
View File

File diff suppressed because it is too large Load Diff