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

pInt not needed anymore

This commit is contained in:
Martin Diehl 2019-05-14 22:44:38 +02:00
parent 387c45d0f4
commit 51e19048f7
15 changed files with 1614 additions and 1717 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
src/DAMASK_marc.f90
View File

@ -96,14 +96,12 @@ end subroutine DAMASK_interface_init
!> @brief solver job name (no extension) as combination of geometry and load case name !> @brief solver job name (no extension) as combination of geometry and load case name
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function getSolverJobName() function getSolverJobName()
use prec, only: & use prec
pReal, &
pInt
implicit none implicit none
character(1024) :: getSolverJobName, inputName character(1024) :: getSolverJobName, inputName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
integer(pInt) :: extPos integer :: extPos
getSolverJobName='' getSolverJobName=''
inputName='' inputName=''
@ -133,9 +131,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
dispt,coord,ffn,frotn,strechn,eigvn,ffn1,frotn1, & dispt,coord,ffn,frotn,strechn,eigvn,ffn1,frotn1, &
strechn1,eigvn1,ncrd,itel,ndeg,ndm,nnode, & strechn1,eigvn1,ncrd,itel,ndeg,ndm,nnode, &
jtype,lclass,ifr,ifu) jtype,lclass,ifr,ifu)
use prec, only: & use prec
pReal, &
pInt
use numerics, only: & use numerics, only: &
!$ DAMASK_NumThreadsInt, & !$ DAMASK_NumThreadsInt, &
numerics_unitlength, & numerics_unitlength, &
@ -180,7 +176,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
implicit none implicit none
!$ include "omp_lib.h" ! the openMP function library !$ include "omp_lib.h" ! the openMP function library
integer(pInt), intent(in) :: & ! according to MSC.Marc 2012 Manual D integer, intent(in) :: & ! according to MSC.Marc 2012 Manual D
ngens, & !< size of stress-strain law ngens, & !< size of stress-strain law
nn, & !< integration point number nn, & !< integration point number
ndi, & !< number of direct components ndi, & !< number of direct components
@ -193,7 +189,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
jtype, & !< element type jtype, & !< element type
ifr, & !< set to 1 if R has been calculated ifr, & !< set to 1 if R has been calculated
ifu !< set to 1 if stretch has been calculated ifu !< set to 1 if stretch has been calculated
integer(pInt), dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D integer, dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D
m, & !< (1) user element number, (2) internal element number m, & !< (1) user element number, (2) internal element number
matus, & !< (1) user material identification number, (2) internal material identification number matus, & !< (1) user material identification number, (2) internal material identification number
kcus, & !< (1) layer number, (2) internal layer number kcus, & !< (1) layer number, (2) internal layer number
@ -236,10 +232,10 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
logical :: cutBack logical :: cutBack
real(pReal), dimension(6) :: stress real(pReal), dimension(6) :: stress
real(pReal), dimension(6,6) :: ddsdde real(pReal), dimension(6,6) :: ddsdde
integer(pInt) :: computationMode, i, cp_en, node, CPnodeID integer :: computationMode, i, cp_en, node, CPnodeID
!$ integer(4) :: defaultNumThreadsInt !< default value set by Marc !$ integer(4) :: defaultNumThreadsInt !< default value set by Marc
if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0_pInt) then if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then
write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn
write(6,'(a,2(i1))') ' Jacobian: ', ngens,ngens write(6,'(a,2(i1))') ' Jacobian: ', ngens,ngens
write(6,'(a,i1)') ' Direct stress: ', ndi write(6,'(a,i1)') ' Direct stress: ', ndi
@ -260,7 +256,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
computationMode = 0_pInt ! save initialization value, since it does not result in any calculation computationMode = 0 ! save initialization value, since it does not result in any calculation
if (lovl == 4 ) then ! jacobian requested by marc if (lovl == 4 ) then ! jacobian requested by marc
if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback
computationMode = CPFEM_RESTOREJACOBIAN computationMode = CPFEM_RESTOREJACOBIAN
@ -307,7 +303,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
if (lastLovl /= lovl) then ! first after ping pong if (lastLovl /= lovl) then ! first after ping pong
call debug_reset() ! resets debugging call debug_reset() ! resets debugging
outdatedFFN1 = .false. outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1_pInt cycleCounter = cycleCounter + 1
mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) ! update cell node coordinates mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) ! update cell node coordinates
call mesh_build_ipCoordinates() ! update ip coordinates call mesh_build_ipCoordinates() ! update ip coordinates
endif endif
@ -324,7 +320,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
lastIncConverged = .false. ! reset flag lastIncConverged = .false. ! reset flag
endif endif
do node = 1,theMesh%elem%nNodes do node = 1,theMesh%elem%nNodes
CPnodeID = mesh_element(4_pInt+node,cp_en) CPnodeID = mesh_element(4+node,cp_en)
mesh_node(1:ndeg,CPnodeID) = mesh_node0(1:ndeg,CPnodeID) + numerics_unitlength * dispt(1:ndeg,node) mesh_node(1:ndeg,CPnodeID) = mesh_node0(1:ndeg,CPnodeID) + numerics_unitlength * dispt(1:ndeg,node)
enddo enddo
endif endif
@ -336,7 +332,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
call debug_info() ! first reports (meaningful) debugging call debug_info() ! first reports (meaningful) debugging
call debug_reset() ! and resets debugging call debug_reset() ! and resets debugging
outdatedFFN1 = .false. outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1_pInt cycleCounter = cycleCounter + 1
mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) ! update cell node coordinates mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) ! update cell node coordinates
call mesh_build_ipCoordinates() ! update ip coordinates call mesh_build_ipCoordinates() ! update ip coordinates
endif endif
@ -376,22 +372,20 @@ end subroutine hypela2
!> @brief calculate internal heat generated due to inelastic energy dissipation !> @brief calculate internal heat generated due to inelastic energy dissipation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine flux(f,ts,n,time) subroutine flux(f,ts,n,time)
use prec, only: & use prec
pReal, &
pInt
use thermal_conduction, only: & use thermal_conduction, only: &
thermal_conduction_getSourceAndItsTangent thermal_conduction_getSourceAndItsTangent
use mesh, only: & use mesh, only: &
mesh_FEasCP mesh_FEasCP
implicit none implicit none
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
ts ts
integer(pInt), dimension(10), intent(in) :: & integer, dimension(10), intent(in) :: &
n n
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
time time
real(pReal), dimension(2), intent(out) :: & real(pReal), dimension(2), intent(out) :: &
f f
call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEasCP('elem',n(1))) call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEasCP('elem',n(1)))
@ -404,9 +398,7 @@ subroutine flux(f,ts,n,time)
!> @details select a variable contour plotting (user subroutine). !> @details select a variable contour plotting (user subroutine).
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine uedinc(inc,incsub) subroutine uedinc(inc,incsub)
use prec, only: & use prec
pReal, &
pInt
use CPFEM, only: & use CPFEM, only: &
CPFEM_results CPFEM_results
@ -424,9 +416,7 @@ end subroutine uedinc
!> @details select a variable contour plotting (user subroutine). !> @details select a variable contour plotting (user subroutine).
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd) subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
use prec, only: & use prec
pReal, &
pInt
use mesh, only: & use mesh, only: &
mesh_FEasCP mesh_FEasCP
use IO, only: & use IO, only: &
@ -436,7 +426,7 @@ subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
materialpoint_sizeResults materialpoint_sizeResults
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
m, & !< element number m, & !< element number
nn, & !< integration point number nn, & !< integration point number
layer, & !< layer number layer, & !< layer number
@ -453,7 +443,7 @@ subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
v !< variable v !< variable
if (jpltcd > materialpoint_sizeResults) call IO_error(700_pInt,jpltcd) ! complain about out of bounds error if (jpltcd > materialpoint_sizeResults) call IO_error(700,jpltcd) ! complain about out of bounds error
v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m)) v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m))
end subroutine plotv end subroutine plotv

52
src/FEsolving.f90
View File

@ -11,8 +11,8 @@ module FEsolving
implicit none implicit none
private private
integer(pInt), public :: & integer, public :: &
restartInc = 1_pInt !< needs description restartInc = 1 !< needs description
logical, public :: & logical, public :: &
symmetricSolver = .false., & !< use a symmetric FEM solver symmetricSolver = .false., & !< use a symmetric FEM solver
@ -20,10 +20,10 @@ module FEsolving
restartRead = .false., & !< restart information to continue calculation from saved state restartRead = .false., & !< restart information to continue calculation from saved state
terminallyIll = .false. !< at least one material point is terminally ill terminallyIll = .false. !< at least one material point is terminally ill
integer(pInt), dimension(:,:), allocatable, public :: & integer, dimension(:,:), allocatable, public :: &
FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP
integer(pInt), dimension(2), public :: & integer, dimension(2), public :: &
FEsolving_execElem !< for ping-pong scheme always whole range, otherwise one specific element FEsolving_execElem !< for ping-pong scheme always whole range, otherwise one specific element
character(len=1024), public :: & character(len=1024), public :: &
@ -61,11 +61,11 @@ subroutine FE_init
implicit none implicit none
#if defined(Marc4DAMASK) || defined(Abaqus) #if defined(Marc4DAMASK) || defined(Abaqus)
integer(pInt), parameter :: & integer, parameter :: &
FILEUNIT = 222_pInt FILEUNIT = 222
integer(pInt) :: j integer :: j
character(len=65536) :: tag, line character(len=65536) :: tag, line
integer(pInt), allocatable, dimension(:) :: chunkPos integer, allocatable, dimension(:) :: chunkPos
#endif #endif
write(6,'(/,a)') ' <<<+- FEsolving init -+>>>' write(6,'(/,a)') ' <<<+- FEsolving init -+>>>'
@ -75,35 +75,35 @@ subroutine FE_init
#if defined(Grid) || defined(FEM) #if defined(Grid) || defined(FEM)
restartInc = interface_RestartInc restartInc = interface_RestartInc
if(restartInc < 0_pInt) then if(restartInc < 0) then
call IO_warning(warning_ID=34_pInt) call IO_warning(warning_ID=34)
restartInc = 0_pInt restartInc = 0
endif endif
restartRead = restartInc > 0_pInt ! only read in if "true" restart requested restartRead = restartInc > 0 ! only read in if "true" restart requested
#else #else
call IO_open_inputFile(FILEUNIT,modelName) call IO_open_inputFile(FILEUNIT,modelName)
rewind(FILEUNIT) rewind(FILEUNIT)
do do
read (FILEUNIT,'(a1024)',END=100) line read (FILEUNIT,'(a1024)',END=100) line
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,chunkPos,1)) ! extract key
select case(tag) select case(tag)
case ('solver') case ('solver')
read (FILEUNIT,'(a1024)',END=100) line ! next line read (FILEUNIT,'(a1024)',END=100) line ! next line
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
symmetricSolver = (IO_intValue(line,chunkPos,2_pInt) /= 1_pInt) symmetricSolver = (IO_intValue(line,chunkPos,2) /= 1)
case ('restart') case ('restart')
read (FILEUNIT,'(a1024)',END=100) line ! next line read (FILEUNIT,'(a1024)',END=100) line ! next line
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
restartWrite = iand(IO_intValue(line,chunkPos,1_pInt),1_pInt) > 0_pInt restartWrite = iand(IO_intValue(line,chunkPos,1),1) > 0
restartRead = iand(IO_intValue(line,chunkPos,1_pInt),2_pInt) > 0_pInt restartRead = iand(IO_intValue(line,chunkPos,1),2) > 0
case ('*restart') case ('*restart')
do j=2_pInt,chunkPos(1) do j=2,chunkPos(1)
restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'write') .or. restartWrite restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'write') .or. restartWrite
restartRead = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'read') .or. restartRead restartRead = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'read') .or. restartRead
enddo enddo
if(restartWrite) then if(restartWrite) then
do j=2_pInt,chunkPos(1) do j=2,chunkPos(1)
restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) /= 'frequency=0') .and. restartWrite restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) /= 'frequency=0') .and. restartWrite
enddo enddo
endif endif
@ -118,11 +118,11 @@ subroutine FE_init
do do
read (FILEUNIT,'(a1024)',END=200) line read (FILEUNIT,'(a1024)',END=200) line
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'restart' & if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'restart' &
.and. IO_lc(IO_stringValue(line,chunkPos,2_pInt)) == 'file' & .and. IO_lc(IO_stringValue(line,chunkPos,2)) == 'file' &
.and. IO_lc(IO_stringValue(line,chunkPos,3_pInt)) == 'job' & .and. IO_lc(IO_stringValue(line,chunkPos,3)) == 'job' &
.and. IO_lc(IO_stringValue(line,chunkPos,4_pInt)) == 'id' ) & .and. IO_lc(IO_stringValue(line,chunkPos,4)) == 'id' ) &
modelName = IO_StringValue(line,chunkPos,6_pInt) modelName = IO_StringValue(line,chunkPos,6)
enddo enddo
#else ! QUESTION: is this meaningful for the spectral/FEM case? #else ! QUESTION: is this meaningful for the spectral/FEM case?
call IO_open_inputFile(FILEUNIT,modelName) call IO_open_inputFile(FILEUNIT,modelName)
@ -130,10 +130,10 @@ subroutine FE_init
do do
read (FILEUNIT,'(a1024)',END=200) line read (FILEUNIT,'(a1024)',END=200) line
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
if (IO_lc(IO_stringValue(line,chunkPos,1_pInt))=='*heading') then if (IO_lc(IO_stringValue(line,chunkPos,1))=='*heading') then
read (FILEUNIT,'(a1024)',END=200) line read (FILEUNIT,'(a1024)',END=200) line
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
modelName = IO_StringValue(line,chunkPos,1_pInt) modelName = IO_StringValue(line,chunkPos,1)
endif endif
enddo enddo
#endif #endif
@ -141,7 +141,7 @@ subroutine FE_init
endif endif
#endif #endif
if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0) then
write(6,'(a21,l1)') ' restart writing: ', restartWrite write(6,'(a21,l1)') ' restart writing: ', restartWrite
write(6,'(a21,l1)') ' restart reading: ', restartRead write(6,'(a21,l1)') ' restart reading: ', restartRead
if (restartRead) write(6,'(a,/)') ' restart Job: '//trim(modelName) if (restartRead) write(6,'(a,/)') ' restart Job: '//trim(modelName)

550
src/IO.f90
View File

File diff suppressed because it is too large Load Diff

40
src/damage_local.f90
View File

@ -9,13 +9,13 @@ module damage_local
implicit none implicit none
private private
integer(pInt), dimension(:,:), allocatable, target, public :: & integer, dimension(:,:), allocatable, target, public :: &
damage_local_sizePostResult !< size of each post result output damage_local_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
damage_local_output !< name of each post result output damage_local_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: & integer, dimension(:), allocatable, target, public :: &
damage_local_Noutput !< number of outputs per instance of this damage damage_local_Noutput !< number of outputs per instance of this damage
enum, bind(c) enum, bind(c)
@ -64,9 +64,9 @@ subroutine damage_local_init
implicit none implicit none
integer(pInt) :: maxNinstance,homog,instance,o,i integer :: maxNinstance,homog,instance,o,i
integer(pInt) :: sizeState integer :: sizeState
integer(pInt) :: NofMyHomog, h integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: & integer(kind(undefined_ID)) :: &
outputID outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
@ -76,13 +76,13 @@ subroutine damage_local_init
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>' write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'
maxNinstance = int(count(damage_type == DAMAGE_local_ID),pInt) maxNinstance = int(count(damage_type == DAMAGE_local_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0) return
allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_local_output (maxval(homogenization_Noutput),maxNinstance)) allocate(damage_local_output (maxval(homogenization_Noutput),maxNinstance))
damage_local_output = '' damage_local_output = ''
allocate(damage_local_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) allocate(damage_local_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(damage_local_Noutput (maxNinstance), source=0_pInt) allocate(damage_local_Noutput (maxNinstance), source=0)
allocate(param(maxNinstance)) allocate(param(maxNinstance))
@ -116,7 +116,7 @@ subroutine damage_local_init
! allocate state arrays ! allocate state arrays
sizeState = 1_pInt sizeState = 1
damageState(homog)%sizeState = sizeState damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_local_sizePostResult(:,instance)) damageState(homog)%sizePostResults = sum(damage_local_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog)) allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
@ -148,14 +148,14 @@ function damage_local_updateState(subdt, ip, el)
damageState damageState
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt subdt
logical, dimension(2) :: & logical, dimension(2) :: &
damage_local_updateState damage_local_updateState
integer(pInt) :: & integer :: &
homog, & homog, &
offset offset
real(pReal) :: & real(pReal) :: &
@ -202,12 +202,12 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
source_damage_anisoductile_getRateAndItsTangent source_damage_anisoductile_getRateAndItsTangent
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer(pInt) :: & integer :: &
phase, & phase, &
grain, & grain, &
source, & source, &
@ -260,26 +260,26 @@ function damage_local_postResults(ip,el)
damage damage
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), dimension(sum(damage_local_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: & real(pReal), dimension(sum(damage_local_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
damage_local_postResults damage_local_postResults
integer(pInt) :: & integer :: &
instance, homog, offset, o, c instance, homog, offset, o, c
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el) offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog) instance = damage_typeInstance(homog)
associate(prm => param(instance)) associate(prm => param(instance))
c = 0_pInt c = 0
outputsLoop: do o = 1_pInt,size(prm%outputID) outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o)) select case(prm%outputID(o))
case (damage_ID) case (damage_ID)
damage_local_postResults(c+1_pInt) = damage(homog)%p(offset) damage_local_postResults(c+1) = damage(homog)%p(offset)
c = c + 1 c = c + 1
end select end select
enddo outputsLoop enddo outputsLoop

56
src/damage_nonlocal.f90
View File

@ -4,19 +4,17 @@
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module damage_nonlocal module damage_nonlocal
use prec, only: & use prec
pReal, &
pInt
implicit none implicit none
private private
integer(pInt), dimension(:,:), allocatable, target, public :: & integer, dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult !< size of each post result output damage_nonlocal_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
damage_nonlocal_output !< name of each post result output damage_nonlocal_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: & integer, dimension(:), allocatable, target, public :: &
damage_nonlocal_Noutput !< number of outputs per instance of this damage damage_nonlocal_Noutput !< number of outputs per instance of this damage
enum, bind(c) enum, bind(c)
@ -64,10 +62,10 @@ subroutine damage_nonlocal_init
implicit none implicit none
integer(pInt) :: maxNinstance,homog,instance,o,i integer :: maxNinstance,homog,instance,o,i
integer(pInt) :: sizeState integer :: sizeState
integer(pInt) :: NofMyHomog, h integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: & integer(kind(undefined_ID)) :: &
outputID outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
@ -75,13 +73,13 @@ subroutine damage_nonlocal_init
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>' write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'
maxNinstance = int(count(damage_type == DAMAGE_nonlocal_ID),pInt) maxNinstance = int(count(damage_type == DAMAGE_nonlocal_ID))
if (maxNinstance == 0_pInt) return if (maxNinstance == 0) return
allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance)) allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance))
damage_nonlocal_output = '' damage_nonlocal_output = ''
allocate(damage_nonlocal_Noutput (maxNinstance), source=0_pInt) allocate(damage_nonlocal_Noutput (maxNinstance), source=0)
allocate(param(maxNinstance)) allocate(param(maxNinstance))
@ -114,7 +112,7 @@ subroutine damage_nonlocal_init
! allocate state arrays ! allocate state arrays
sizeState = 1_pInt sizeState = 1
damageState(homog)%sizeState = sizeState damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance)) damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog)) allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
@ -155,12 +153,12 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
source_damage_anisoductile_getRateAndItsTangent source_damage_anisoductile_getRateAndItsTangent
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer(pInt) :: & integer :: &
phase, & phase, &
grain, & grain, &
source, & source, &
@ -218,12 +216,12 @@ function damage_nonlocal_getDiffusion33(ip,el)
crystallite_push33ToRef crystallite_push33ToRef
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
damage_nonlocal_getDiffusion33 damage_nonlocal_getDiffusion33
integer(pInt) :: & integer :: &
homog, & homog, &
grain grain
@ -235,7 +233,7 @@ function damage_nonlocal_getDiffusion33(ip,el)
enddo enddo
damage_nonlocal_getDiffusion33 = & damage_nonlocal_getDiffusion33 = &
charLength**2_pInt*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal) charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
end function damage_nonlocal_getDiffusion33 end function damage_nonlocal_getDiffusion33
@ -252,10 +250,10 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
homogenization_Ngrains homogenization_Ngrains
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
integer(pInt) :: & integer :: &
ipc ipc
damage_nonlocal_getMobility = 0.0_pReal damage_nonlocal_getMobility = 0.0_pReal
@ -279,12 +277,12 @@ subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
damage damage
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer(pInt) :: & integer :: &
homog, & homog, &
offset offset
@ -305,26 +303,26 @@ function damage_nonlocal_postResults(ip,el)
damage damage
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: & real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
damage_nonlocal_postResults damage_nonlocal_postResults
integer(pInt) :: & integer :: &
instance, homog, offset, o, c instance, homog, offset, o, c
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el) offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog) instance = damage_typeInstance(homog)
associate(prm => param(instance)) associate(prm => param(instance))
c = 0_pInt c = 0
outputsLoop: do o = 1_pInt,size(prm%outputID) outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o)) select case(prm%outputID(o))
case (damage_ID) case (damage_ID)
damage_nonlocal_postResults(c+1_pInt) = damage(homog)%p(offset) damage_nonlocal_postResults(c+1) = damage(homog)%p(offset)
c = c + 1 c = c + 1
end select end select
enddo outputsLoop enddo outputsLoop

276
src/grid/DAMASK_grid.f90
View File

@ -9,69 +9,19 @@
program DAMASK_spectral program DAMASK_spectral
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
use PETScsys use PETScsys
use prec, only: & use prec
pInt, & use DAMASK_interface
pLongInt, & use IO
pReal, & use config
tol_math_check, & use debug
dNeq use math
use DAMASK_interface, only: & use mesh
DAMASK_interface_init, & use CPFEM2
loadCaseFile, & use FEsolving
geometryFile, & use numerics
getSolverJobName, & use homogenization
interface_restartInc use material
use IO, only: & use spectral_utilities
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_error, &
IO_lc, &
IO_intOut, &
IO_warning
use config, only: &
config_numerics
use debug, only: &
debug_level, &
debug_spectral, &
debug_levelBasic
use math ! need to include the whole module for FFTW
use mesh, only: &
grid, &
geomSize
use CPFEM2, only: &
CPFEM_initAll, &
CPFEM_results
use FEsolving, only: &
restartWrite, &
restartInc
use numerics, only: &
worldrank, &
worldsize, &
stagItMax, &
maxCutBack, &
continueCalculation
use homogenization, only: &
materialpoint_sizeResults, &
materialpoint_results, &
materialpoint_postResults
use material, only: &
thermal_type, &
damage_type, &
THERMAL_conduction_ID, &
DAMAGE_nonlocal_ID
use spectral_utilities, only: &
utilities_init, &
tSolutionState, &
tLoadCase, &
cutBack, &
nActiveFields, &
FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID
use grid_mech_spectral_basic use grid_mech_spectral_basic
use grid_mech_spectral_polarisation use grid_mech_spectral_polarisation
use grid_mech_FEM use grid_mech_FEM
@ -86,11 +36,11 @@ program DAMASK_spectral
! variables related to information from load case and geom file ! variables related to information from load case and geom file
real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0) real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors
integer(pInt), allocatable, dimension(:) :: chunkPos integer, allocatable, dimension(:) :: chunkPos
integer(pInt) :: & integer :: &
N_t = 0_pInt, & !< # of time indicators found in load case file N_t = 0, & !< # of time indicators found in load case file
N_n = 0_pInt, & !< # of increment specifiers found in load case file N_n = 0, & !< # of increment specifiers found in load case file
N_def = 0_pInt !< # of rate of deformation specifiers found in load case file N_def = 0 !< # of rate of deformation specifiers found in load case file
character(len=65536) :: & character(len=65536) :: &
line line
@ -99,8 +49,8 @@ program DAMASK_spectral
real(pReal), dimension(3,3), parameter :: & real(pReal), dimension(3,3), parameter :: &
ones = 1.0_pReal, & ones = 1.0_pReal, &
zeros = 0.0_pReal zeros = 0.0_pReal
integer(pInt), parameter :: & integer, parameter :: &
subStepFactor = 2_pInt !< for each substep, divide the last time increment by 2.0 subStepFactor = 2 !< for each substep, divide the last time increment by 2.0
real(pReal) :: & real(pReal) :: &
time = 0.0_pReal, & !< elapsed time time = 0.0_pReal, & !< elapsed time
time0 = 0.0_pReal, & !< begin of interval time0 = 0.0_pReal, & !< begin of interval
@ -110,21 +60,21 @@ program DAMASK_spectral
logical :: & logical :: &
guess, & !< guess along former trajectory guess, & !< guess along former trajectory
stagIterate stagIterate
integer(pInt) :: & integer :: &
i, j, k, l, field, & i, j, k, l, field, &
errorID = 0_pInt, & errorID = 0, &
cutBackLevel = 0_pInt, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0_pInt !< fraction of current time interval stepFraction = 0 !< fraction of current time interval
integer(pInt) :: & integer :: &
currentLoadcase = 0_pInt, & !< current load case currentLoadcase = 0, & !< current load case
inc, & !< current increment in current load case inc, & !< current increment in current load case
totalIncsCounter = 0_pInt, & !< total # of increments totalIncsCounter = 0, & !< total # of increments
convergedCounter = 0_pInt, & !< # of converged increments convergedCounter = 0, & !< # of converged increments
notConvergedCounter = 0_pInt, & !< # of non-converged increments notConvergedCounter = 0, & !< # of non-converged increments
fileUnit = 0_pInt, & !< file unit for reading load case and writing results fileUnit = 0, & !< file unit for reading load case and writing results
myStat, & myStat, &
statUnit = 0_pInt, & !< file unit for statistics output statUnit = 0, & !< file unit for statistics output
lastRestartWritten = 0_pInt, & !< total increment # at which last restart information was written lastRestartWritten = 0, & !< total increment # at which last restart information was written
stagIter stagIter
character(len=6) :: loadcase_string character(len=6) :: loadcase_string
character(len=1024) :: & character(len=1024) :: &
@ -134,8 +84,8 @@ program DAMASK_spectral
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
integer(MPI_OFFSET_KIND) :: fileOffset integer(MPI_OFFSET_KIND) :: fileOffset
integer(MPI_OFFSET_KIND), dimension(:), allocatable :: outputSize integer(MPI_OFFSET_KIND), dimension(:), allocatable :: outputSize
integer(pInt), parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742 integer, parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742
integer(pInt), parameter :: maxRealOut = maxByteOut/pReal integer, parameter :: maxRealOut = maxByteOut/pReal
integer(pLongInt), dimension(2) :: outputIndex integer(pLongInt), dimension(2) :: outputIndex
PetscErrorCode :: ierr PetscErrorCode :: ierr
procedure(grid_mech_spectral_basic_init), pointer :: & procedure(grid_mech_spectral_basic_init), pointer :: &
@ -174,20 +124,20 @@ program DAMASK_spectral
case ('polarisation') case ('polarisation')
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence') call IO_warning(42, ext_msg='debug Divergence')
mech_init => grid_mech_spectral_polarisation_init mech_init => grid_mech_spectral_polarisation_init
mech_forward => grid_mech_spectral_polarisation_forward mech_forward => grid_mech_spectral_polarisation_forward
mech_solution => grid_mech_spectral_polarisation_solution mech_solution => grid_mech_spectral_polarisation_solution
case ('fem') case ('fem')
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence') call IO_warning(42, ext_msg='debug Divergence')
mech_init => grid_mech_FEM_init mech_init => grid_mech_FEM_init
mech_forward => grid_mech_FEM_forward mech_forward => grid_mech_FEM_forward
mech_solution => grid_mech_FEM_solution mech_solution => grid_mech_FEM_solution
case default case default
call IO_error(error_ID = 891_pInt, ext_msg = config_numerics%getString('spectral_solver')) call IO_error(error_ID = 891, ext_msg = config_numerics%getString('spectral_solver'))
end select end select
@ -195,27 +145,27 @@ program DAMASK_spectral
! reading information from load case file and to sanity checks ! reading information from load case file and to sanity checks
allocate (loadCases(0)) ! array of load cases allocate (loadCases(0)) ! array of load cases
open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read') open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read')
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=trim(loadCaseFile)) if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=trim(loadCaseFile))
do do
read(fileUnit, '(A)', iostat=myStat) line read(fileUnit, '(A)', iostat=myStat) line
if ( myStat /= 0_pInt) exit if ( myStat /= 0) exit
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
currentLoadCase = currentLoadCase + 1_pInt currentLoadCase = currentLoadCase + 1
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1) ! reading compulsory parameters for loadcase do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i))) select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('l','velocitygrad','velgrad','velocitygradient','fdot','dotf','f') case('l','velocitygrad','velgrad','velocitygradient','fdot','dotf','f')
N_def = N_def + 1_pInt N_def = N_def + 1
case('t','time','delta') case('t','time','delta')
N_t = N_t + 1_pInt N_t = N_t + 1
case('n','incs','increments','steps','logincs','logincrements','logsteps') case('n','incs','increments','steps','logincs','logincrements','logsteps')
N_n = N_n + 1_pInt N_n = N_n + 1
end select end select
enddo enddo
if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1_pInt) & ! sanity check if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check
call IO_error(error_ID=837_pInt,el=currentLoadCase,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase
newLoadCase%stress%myType='stress' newLoadCase%stress%myType='stress'
field = 1 field = 1
@ -229,7 +179,7 @@ program DAMASK_spectral
newLoadCase%ID(field) = FIELD_DAMAGE_ID newLoadCase%ID(field) = FIELD_DAMAGE_ID
endif damageActive endif damageActive
readIn: do i = 1_pInt, chunkPos(1) readIn: do i = 1, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i))) select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient','f') ! assign values for the deformation BC matrix case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient','f') ! assign values for the deformation BC matrix
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
@ -241,7 +191,7 @@ program DAMASK_spectral
else else
newLoadCase%deformation%myType = 'l' newLoadCase%deformation%myType = 'l'
endif endif
do j = 1_pInt, 9_pInt do j = 1, 9
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a * temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a *
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo enddo
@ -250,7 +200,7 @@ program DAMASK_spectral
newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation
case('p','pk1','piolakirchhoff','stress', 's') case('p','pk1','piolakirchhoff','stress', 's')
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
do j = 1_pInt, 9_pInt do j = 1, 9
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo enddo
@ -258,54 +208,54 @@ program DAMASK_spectral
newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical) newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical)
newLoadCase%stress%values = math_9to33(temp_valueVector) newLoadCase%stress%values = math_9to33(temp_valueVector)
case('t','time','delta') ! increment time case('t','time','delta') ! increment time
newLoadCase%time = IO_floatValue(line,chunkPos,i+1_pInt) newLoadCase%time = IO_floatValue(line,chunkPos,i+1)
case('n','incs','increments','steps') ! number of increments case('n','incs','increments','steps') ! number of increments
newLoadCase%incs = IO_intValue(line,chunkPos,i+1_pInt) newLoadCase%incs = IO_intValue(line,chunkPos,i+1)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
newLoadCase%incs = IO_intValue(line,chunkPos,i+1_pInt) newLoadCase%incs = IO_intValue(line,chunkPos,i+1)
newLoadCase%logscale = 1_pInt newLoadCase%logscale = 1
case('freq','frequency','outputfreq') ! frequency of result writings case('freq','frequency','outputfreq') ! frequency of result writings
newLoadCase%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt) newLoadCase%outputfrequency = IO_intValue(line,chunkPos,i+1)
case('r','restart','restartwrite') ! frequency of writing restart information case('r','restart','restartwrite') ! frequency of writing restart information
newLoadCase%restartfrequency = & newLoadCase%restartfrequency = &
max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt)) max(0,IO_intValue(line,chunkPos,i+1))
case('guessreset','dropguessing') case('guessreset','dropguessing')
newLoadCase%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory newLoadCase%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('euler') ! rotation of load case given in euler angles case('euler') ! rotation of load case given in euler angles
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
l = 1_pInt ! assuming values given in degrees l = 1 ! assuming values given in degrees
k = 1_pInt ! assuming keyword indicating degree/radians present k = 1 ! assuming keyword indicating degree/radians present
select case (IO_lc(IO_stringValue(line,chunkPos,i+1_pInt))) select case (IO_lc(IO_stringValue(line,chunkPos,i+1)))
case('deg','degree') case('deg','degree')
case('rad','radian') ! don't convert from degree to radian case('rad','radian') ! don't convert from degree to radian
l = 0_pInt l = 0
case default case default
k = 0_pInt k = 0
end select end select
do j = 1_pInt, 3_pInt do j = 1, 3
temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j)
enddo enddo
if (l == 1_pInt) temp_valueVector(1:3) = temp_valueVector(1:3) * inRad ! convert to rad if (l == 1) temp_valueVector(1:3) = temp_valueVector(1:3) * INRAD ! convert to rad
newLoadCase%rotation = math_EulerToR(temp_valueVector(1:3)) ! convert rad Eulers to rotation matrix newLoadCase%rotation = math_EulerToR(temp_valueVector(1:3)) ! convert rad Eulers to rotation matrix
case('rotation','rot') ! assign values for the rotation matrix case('rotation','rot') ! assign values for the rotation matrix
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
do j = 1_pInt, 9_pInt do j = 1, 9
temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j)
enddo enddo
newLoadCase%rotation = math_9to33(temp_valueVector) newLoadCase%rotation = math_9to33(temp_valueVector)
end select end select
enddo readIn enddo readIn
newLoadCase%followFormerTrajectory = merge(.true.,.false.,currentLoadCase > 1_pInt) ! by default, guess from previous load case newLoadCase%followFormerTrajectory = merge(.true.,.false.,currentLoadCase > 1) ! by default, guess from previous load case
reportAndCheck: if (worldrank == 0) then reportAndCheck: if (worldrank == 0) then
write (loadcase_string, '(i6)' ) currentLoadCase write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(/,1x,a,i6)') 'load case: ', currentLoadCase write(6,'(/,1x,a,i6)') 'load case: ', currentLoadCase
if (.not. newLoadCase%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory' if (.not. newLoadCase%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory'
if (newLoadCase%deformation%myType == 'l') then if (newLoadCase%deformation%myType == 'l') then
do j = 1_pInt, 3_pInt do j = 1, 3
if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. &
any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined
enddo enddo
write(6,'(2x,a)') 'velocity gradient:' write(6,'(2x,a)') 'velocity gradient:'
else if (newLoadCase%deformation%myType == 'f') then else if (newLoadCase%deformation%myType == 'f') then
@ -313,7 +263,7 @@ program DAMASK_spectral
else else
write(6,'(2x,a)') 'deformation gradient rate:' write(6,'(2x,a)') 'deformation gradient rate:'
endif endif
do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt do i = 1, 3; do j = 1, 3
if(newLoadCase%deformation%maskLogical(i,j)) then if(newLoadCase%deformation%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j) write(6,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j)
else else
@ -322,13 +272,13 @@ program DAMASK_spectral
enddo; write(6,'(/)',advance='no') enddo; write(6,'(/)',advance='no')
enddo enddo
if (any(newLoadCase%stress%maskLogical .eqv. & if (any(newLoadCase%stress%maskLogical .eqv. &
newLoadCase%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only
if (any(newLoadCase%stress%maskLogical .and. & if (any(newLoadCase%stress%maskLogical .and. &
transpose(newLoadCase%stress%maskLogical) .and. & transpose(newLoadCase%stress%maskLogical) .and. &
reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
errorID = 838_pInt ! no rotation is allowed by stress BC errorID = 838 ! no rotation is allowed by stress BC
write(6,'(2x,a)') 'stress / GPa:' write(6,'(2x,a)') 'stress / GPa:'
do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt do i = 1, 3; do j = 1, 3
if(newLoadCase%stress%maskLogical(i,j)) then if(newLoadCase%stress%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal write(6,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal
else else
@ -340,18 +290,18 @@ program DAMASK_spectral
transpose(newLoadCase%rotation))-math_I3) > & transpose(newLoadCase%rotation))-math_I3) > &
reshape(spread(tol_math_check,1,9),[ 3,3]))& reshape(spread(tol_math_check,1,9),[ 3,3]))&
.or. abs(math_det33(newLoadCase%rotation)) > & .or. abs(math_det33(newLoadCase%rotation)) > &
1.0_pReal + tol_math_check) errorID = 846_pInt ! given rotation matrix contains strain 1.0_pReal + tol_math_check) errorID = 846 ! given rotation matrix contains strain
if (any(dNeq(newLoadCase%rotation, math_I3))) & if (any(dNeq(newLoadCase%rotation, math_I3))) &
write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
transpose(newLoadCase%rotation) transpose(newLoadCase%rotation)
if (newLoadCase%time < 0.0_pReal) errorID = 834_pInt ! negative time increment if (newLoadCase%time < 0.0_pReal) errorID = 834 ! negative time increment
write(6,'(2x,a,f12.6)') 'time: ', newLoadCase%time write(6,'(2x,a,f12.6)') 'time: ', newLoadCase%time
if (newLoadCase%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count if (newLoadCase%incs < 1) errorID = 835 ! non-positive incs count
write(6,'(2x,a,i5)') 'increments: ', newLoadCase%incs write(6,'(2x,a,i5)') 'increments: ', newLoadCase%incs
if (newLoadCase%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency if (newLoadCase%outputfrequency < 1) errorID = 836 ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', newLoadCase%outputfrequency write(6,'(2x,a,i5)') 'output frequency: ', newLoadCase%outputfrequency
write(6,'(2x,a,i5)') 'restart frequency: ', newLoadCase%restartfrequency write(6,'(2x,a,i5)') 'restart frequency: ', newLoadCase%restartfrequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
endif reportAndCheck endif reportAndCheck
loadCases = [loadCases,newLoadCase] ! load case is ok, append it loadCases = [loadCases,newLoadCase] ! load case is ok, append it
enddo enddo
@ -383,7 +333,7 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write header of output file ! write header of output file
if (worldrank == 0) then if (worldrank == 0) then
writeHeader: if (interface_restartInc < 1_pInt) then writeHeader: if (interface_restartInc < 1) then
open(newunit=fileUnit,file=trim(getSolverJobName())//& open(newunit=fileUnit,file=trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED',status='REPLACE') '.spectralOut',form='UNFORMATTED',status='REPLACE')
write(fileUnit) 'load:', trim(loadCaseFile) ! ... and write header write(fileUnit) 'load:', trim(loadCaseFile) ! ... and write header
@ -417,59 +367,59 @@ program DAMASK_spectral
allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND) allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND)
outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND) outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND)
call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_allreduce') if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_allreduce')
call MPI_file_open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.spectralOut', & call MPI_file_open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.spectralOut', &
MPI_MODE_WRONLY + MPI_MODE_APPEND, & MPI_MODE_WRONLY + MPI_MODE_APPEND, &
MPI_INFO_NULL, & MPI_INFO_NULL, &
fileUnit, & fileUnit, &
ierr) ierr)
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_open') if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_open')
call MPI_file_get_position(fileUnit,fileOffset,ierr) ! get offset from header call MPI_file_get_position(fileUnit,fileOffset,ierr) ! get offset from header
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_get_position') if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_get_position')
fileOffset = fileOffset + sum(outputSize(1:worldrank)) ! offset of my process in file (header + processes before me) fileOffset = fileOffset + sum(outputSize(1:worldrank)) ! offset of my process in file (header + processes before me)
call MPI_file_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr) call MPI_file_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr)
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_seek') if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_seek')
writeUndeformed: if (interface_restartInc < 1_pInt) then writeUndeformed: if (interface_restartInc < 1) then
write(6,'(1/,a)') ' ... writing initial configuration to file ........................' write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
call CPFEM_results(0_pInt,0.0_pReal) call CPFEM_results(0,0.0_pReal)
do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & ! QUESTION: why not starting i at 0 instead of murky 1? outputIndex = int([(i-1)*((maxRealOut)/materialpoint_sizeResults)+1, &
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), & call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
[(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), & [(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)), & int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)), &
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr) MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write') if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_write')
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif writeUndeformed endif writeUndeformed
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) loadCaseLooping: do currentLoadCase = 1, size(loadCases)
time0 = time ! load case start time time0 = time ! load case start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs incLooping: do inc = 1, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt totalIncsCounter = totalIncsCounter + 1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forwarding time ! forwarding time
timeIncOld = timeinc ! last timeinc that brought former inc to an end timeIncOld = timeinc ! last timeinc that brought former inc to an end
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale if (loadCases(currentLoadCase)%logscale == 0) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else else
if (currentLoadCase == 1_pInt) then ! 1st load case of logarithmic scale if (currentLoadCase == 1) then ! 1st load case of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st load case of logarithmic scale if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st load case of logarithmic scale else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal)) timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal))
endif endif
else ! not-1st load case of logarithmic scale else ! not-1st load case of logarithmic scale
timeinc = time0 * & timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/& ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))& real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1_pInt ,pReal)/& -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1 ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))) real(loadCases(currentLoadCase)%incs ,pReal)))
endif endif
endif endif
@ -479,12 +429,12 @@ program DAMASK_spectral
time = time + timeinc ! just advance time, skip already performed calculation time = time + timeinc ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping else skipping
stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time time = time + timeinc ! forward target time
stepFraction = stepFraction + 1_pInt ! count step stepFraction = stepFraction + 1 ! count step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
@ -524,7 +474,7 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve fields ! solve fields
stagIter = 0_pInt stagIter = 0
stagIterate = .true. stagIterate = .true.
do while (stagIterate) do while (stagIterate)
do field = 1, nActiveFields do field = 1, nActiveFields
@ -546,7 +496,7 @@ program DAMASK_spectral
if (.not. solres(field)%converged) exit ! no solution found if (.not. solres(field)%converged) exit ! no solution found
enddo enddo
stagIter = stagIter + 1_pInt stagIter = stagIter + 1
stagIterate = stagIter < stagItMax & stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) & .and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
@ -567,52 +517,52 @@ program DAMASK_spectral
endif endif
elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated? elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated?
cutBack = .true. cutBack = .true.
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time time = time - timeinc ! rewind time
timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep
write(6,'(/,a)') ' cutting back ' write(6,'(/,a)') ' cutting back '
else ! no more options to continue else ! no more options to continue
call IO_warning(850_pInt) call IO_warning(850)
call MPI_file_close(fileUnit,ierr) call MPI_file_close(fileUnit,ierr)
close(statUnit) close(statUnit)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written call quit(-1*(lastRestartWritten+1)) ! quit and provide information about last restart inc written
endif endif
enddo subStepLooping enddo subStepLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then if (all(solres(:)%converged)) then
convergedCounter = convergedCounter + 1_pInt convergedCounter = convergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
' increment ', totalIncsCounter, ' converged' ' increment ', totalIncsCounter, ' converged'
else else
notConvergedCounter = notConvergedCounter + 1_pInt notConvergedCounter = notConvergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged' ' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6) endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................' write(6,'(1/,a)') ' ... writing results to file ......................................'
flush(6) flush(6)
call materialpoint_postResults() call materialpoint_postResults()
call MPI_file_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr) call MPI_file_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr)
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek') if (ierr /= 0) call IO_error(894, ext_msg='MPI_file_seek')
do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & outputIndex=int([(i-1)*((maxRealOut)/materialpoint_sizeResults)+1, &
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),& call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),&
[(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), & [(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)),& int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)),&
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr) MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write') if(ierr /=0) call IO_error(894, ext_msg='MPI_file_write')
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
call CPFEM_results(totalIncsCounter,time) call CPFEM_results(totalIncsCounter,time)
endif endif
if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ... if ( loadCases(currentLoadCase)%restartFrequency > 0 & ! writing of restart info requested ...
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information .and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then ! ... and at frequency of writing restart information
restartWrite = .true. ! set restart parameter for FEsolving restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write? lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif endif
@ -636,7 +586,7 @@ program DAMASK_spectral
call MPI_file_close(fileUnit,ierr) call MPI_file_close(fileUnit,ierr)
close(statUnit) close(statUnit)
if (notConvergedCounter > 0_pInt) call quit(2_pInt) ! error if some are not converged if (notConvergedCounter > 0) call quit(2) ! error if some are not converged
call quit(0_pInt) ! no complains ;) call quit(0) ! no complains ;)
end program DAMASK_spectral end program DAMASK_spectral

68
src/kinematics_cleavage_opening.f90
View File

@ -5,18 +5,16 @@
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module kinematics_cleavage_opening module kinematics_cleavage_opening
use prec, only: & use prec
pReal, &
pInt
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, private :: kinematics_cleavage_opening_instance integer, dimension(:), allocatable, private :: kinematics_cleavage_opening_instance
type, private :: tParameters !< container type for internal constitutive parameters type, private :: tParameters !< container type for internal constitutive parameters
integer(pInt) :: & integer :: &
totalNcleavage totalNcleavage
integer(pInt), dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
Ncleavage !< active number of cleavage systems per family Ncleavage !< active number of cleavage systems per family
real(pReal) :: & real(pReal) :: &
sdot0, & sdot0, &
@ -27,17 +25,17 @@ module kinematics_cleavage_opening
end type end type
! Begin Deprecated ! Begin Deprecated
integer(pInt), dimension(:), allocatable, private :: & integer, dimension(:), allocatable, private :: &
kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems
integer(pInt), dimension(:,:), allocatable, private :: & integer, dimension(:,:), allocatable, private :: &
kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family
real(pReal), dimension(:), allocatable, private :: & real(pReal), dimension(:), allocatable, private :: &
kinematics_cleavage_opening_sdot_0, & kinematics_cleavage_opening_sdot_0, &
kinematics_cleavage_opening_N kinematics_cleavage_opening_N
real(pReal), dimension(:,:), allocatable, private :: & real(pReal), dimension(:,:), allocatable, private :: &
kinematics_cleavage_opening_critDisp, & kinematics_cleavage_opening_critDisp, &
kinematics_cleavage_opening_critLoad kinematics_cleavage_opening_critLoad
! End Deprecated ! End Deprecated
@ -71,32 +69,32 @@ subroutine kinematics_cleavage_opening_init()
lattice_NcleavageSystem lattice_NcleavageSystem
implicit none implicit none
integer(pInt), allocatable, dimension(:) :: tempInt integer, allocatable, dimension(:) :: tempInt
real(pReal), allocatable, dimension(:) :: tempFloat real(pReal), allocatable, dimension(:) :: tempFloat
integer(pInt) :: maxNinstance,p,instance,kinematics integer :: maxNinstance,p,instance,kinematics
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>' write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'
maxNinstance = int(count(phase_kinematics == KINEMATICS_cleavage_opening_ID),pInt) maxNinstance = int(count(phase_kinematics == KINEMATICS_cleavage_opening_ID))
if (maxNinstance == 0_pInt) return if (maxNinstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0_pInt) allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0)
do p = 1_pInt, size(config_phase) do p = 1, size(config_phase)
kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct? kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct?
enddo enddo
allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt) allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0)
allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0_pInt) allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0)
allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal) allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal) allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal)
do p = 1_pInt, size(config_phase) do p = 1, size(config_phase)
if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle
instance = kinematics_cleavage_opening_instance(p) instance = kinematics_cleavage_opening_instance(p)
kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0') kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0')
@ -115,13 +113,13 @@ subroutine kinematics_cleavage_opening_init()
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance)) kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance))
kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether
if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) & if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')') call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) & if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')') call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) & if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')') call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) & if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')') call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
enddo enddo
end subroutine kinematics_cleavage_opening_init end subroutine kinematics_cleavage_opening_init
@ -130,8 +128,6 @@ end subroutine kinematics_cleavage_opening_init
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
use prec, only: &
tol_math_check
use math, only: & use math, only: &
math_mul33xx33 math_mul33xx33
use material, only: & use material, only: &
@ -145,7 +141,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
lattice_NcleavageSystem lattice_NcleavageSystem
implicit none implicit none
integer(pInt), intent(in) :: & integer, intent(in) :: &
ipc, & !< grain number ipc, & !< grain number
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
@ -155,7 +151,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld !< damage velocity gradient Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
integer(pInt) :: & integer :: &
instance, phase, & instance, phase, &
homog, damageOffset, & homog, damageOffset, &
f, i, index_myFamily, k, l, m, n f, i, index_myFamily, k, l, m, n
@ -170,9 +166,9 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld = 0.0_pReal Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal dLd_dTstar = 0.0_pReal
do f = 1_pInt,lattice_maxNcleavageFamily do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1_pInt,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family do i = 1,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
@ -186,7 +182,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase) Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)
dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ & dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_d) - traction_crit) max(0.0_pReal, abs(traction_d) - traction_crit)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* & dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,1,index_myFamily+i,phase) lattice_Scleavage(m,n,1,index_myFamily+i,phase)
@ -200,7 +196,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase) Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)
dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ & dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_t) - traction_crit) max(0.0_pReal, abs(traction_t) - traction_crit)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* & dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,2,index_myFamily+i,phase) lattice_Scleavage(m,n,2,index_myFamily+i,phase)
@ -214,7 +210,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase) Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ & dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_n) - traction_crit) max(0.0_pReal, abs(traction_n) - traction_crit)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* & dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,3,index_myFamily+i,phase) lattice_Scleavage(m,n,3,index_myFamily+i,phase)

42
src/kinematics_slipplane_opening.f90
View File

@ -5,18 +5,16 @@
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module kinematics_slipplane_opening module kinematics_slipplane_opening
use prec, only: & use prec
pReal, &
pInt
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, private :: kinematics_slipplane_opening_instance integer, dimension(:), allocatable, private :: kinematics_slipplane_opening_instance
type, private :: tParameters !< container type for internal constitutive parameters type, private :: tParameters !< container type for internal constitutive parameters
integer(pInt) :: & integer :: &
totalNslip totalNslip
integer(pInt), dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
Nslip !< active number of slip systems per family Nslip !< active number of slip systems per family
real(pReal) :: & real(pReal) :: &
sdot0, & sdot0, &
@ -58,26 +56,25 @@ subroutine kinematics_slipplane_opening_init()
KINEMATICS_slipplane_opening_ID KINEMATICS_slipplane_opening_ID
use lattice use lattice
implicit none
integer(pInt) :: maxNinstance,p,instance,kinematics integer :: maxNinstance,p,instance,kinematics
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>' write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'
maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID) maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID)
if (maxNinstance == 0) return if (maxNinstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0_pInt) allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0)
do p = 1_pInt, size(config_phase) do p = 1, size(config_phase)
kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct? kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct?
enddo enddo
allocate(param(maxNinstance)) allocate(param(maxNinstance))
do p = 1_pInt, size(config_phase) do p = 1, size(config_phase)
if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle
associate(prm => param(kinematics_slipplane_opening_instance(p)), & associate(prm => param(kinematics_slipplane_opening_instance(p)), &
config => config_phase(p)) config => config_phase(p))
@ -91,19 +88,19 @@ subroutine kinematics_slipplane_opening_init()
prm%critLoad = math_expand(prm%critLoad, prm%Nslip) prm%critLoad = math_expand(prm%critLoad, prm%Nslip)
prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),& prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),& prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),& prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
! if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) & ! if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) &
! call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')') ! call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) & ! if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) &
! call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')') ! call IO_error(211,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) & ! if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) &
! call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')') ! call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')')
end associate end associate
enddo enddo
@ -114,8 +111,6 @@ end subroutine kinematics_slipplane_opening_init
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
use prec, only: &
tol_math_check
use math, only: & use math, only: &
math_mul33xx33, & math_mul33xx33, &
math_outer math_outer
@ -125,7 +120,6 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
damage, & damage, &
damageMapping damageMapping
implicit none
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< grain number ipc, & !< grain number
ip, & !< integration point number ip, & !< integration point number
@ -173,7 +167,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
if (abs(udotd) > tol_math_check) then if (abs(udotd) > tol_math_check) then
Ld = Ld + udotd*projection_d Ld = Ld + udotd*projection_d
dudotd_dt = udotd*prm%n/traction_d dudotd_dt = udotd*prm%n/traction_d
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotd_dt*projection_d(k,l)*projection_d(m,n) dudotd_dt*projection_d(k,l)*projection_d(m,n)
endif endif
@ -185,7 +179,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
if (abs(udott) > tol_math_check) then if (abs(udott) > tol_math_check) then
Ld = Ld + udott*projection_t Ld = Ld + udott*projection_t
dudott_dt = udott*prm%n/traction_t dudott_dt = udott*prm%n/traction_t
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudott_dt*projection_t(k,l)*projection_t(m,n) dudott_dt*projection_t(k,l)*projection_t(m,n)
endif endif
@ -197,7 +191,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
if (abs(udotn) > tol_math_check) then if (abs(udotn) > tol_math_check) then
Ld = Ld + udotn*projection_n Ld = Ld + udotn*projection_n
dudotn_dt = udotn*prm%n/traction_n dudotn_dt = udotn*prm%n/traction_n
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotn_dt*projection_n(k,l)*projection_n(m,n) dudotn_dt*projection_n(k,l)*projection_n(m,n)
endif endif

234
src/material.f90
View File

@ -98,10 +98,10 @@ module material
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: &
homogenization_type !< type of each homogenization homogenization_type !< type of each homogenization
integer(pInt), public, protected :: & integer, public, protected :: &
homogenization_maxNgrains !< max number of grains in any USED homogenization homogenization_maxNgrains !< max number of grains in any USED homogenization
integer(pInt), dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
phase_Nsources, & !< number of source mechanisms active in each phase phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
@ -132,7 +132,7 @@ module material
! END NEW MAPPINGS ! END NEW MAPPINGS
! DEPRECATED: use material_phaseAt ! DEPRECATED: use material_phaseAt
integer(pInt), dimension(:,:,:), allocatable, public :: & integer, dimension(:,:,:), allocatable, public :: &
material_phase !< phase (index) of each grain,IP,element material_phase !< phase (index) of each grain,IP,element
type(tPlasticState), allocatable, dimension(:), public :: & type(tPlasticState), allocatable, dimension(:), public :: &
@ -144,7 +144,7 @@ module material
thermalState, & thermalState, &
damageState damageState
integer(pInt), dimension(:,:,:), allocatable, public, protected :: & integer, dimension(:,:,:), allocatable, public, protected :: &
material_texture !< texture (index) of each grain,IP,element material_texture !< texture (index) of each grain,IP,element
real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:), allocatable, public, protected :: &
@ -155,15 +155,15 @@ module material
microstructure_elemhomo, & !< flag to indicate homogeneous microstructure distribution over element's IPs microstructure_elemhomo, & !< flag to indicate homogeneous microstructure distribution over element's IPs
phase_localPlasticity !< flags phases with local constitutive law phase_localPlasticity !< flags phases with local constitutive law
integer(pInt), private :: & integer, private :: &
microstructure_maxNconstituents, & !< max number of constituents in any phase microstructure_maxNconstituents, & !< max number of constituents in any phase
texture_maxNgauss !< max number of Gauss components in any texture texture_maxNgauss !< max number of Gauss components in any texture
integer(pInt), dimension(:), allocatable, private :: & integer, dimension(:), allocatable, private :: &
microstructure_Nconstituents, & !< number of constituents in each microstructure microstructure_Nconstituents, & !< number of constituents in each microstructure
texture_Ngauss !< number of Gauss components per texture texture_Ngauss !< number of Gauss components per texture
integer(pInt), dimension(:,:), allocatable, private :: & integer, dimension(:,:), allocatable, private :: &
microstructure_phase, & !< phase IDs of each microstructure microstructure_phase, & !< phase IDs of each microstructure
microstructure_texture !< texture IDs of each microstructure microstructure_texture !< texture IDs of each microstructure
@ -178,11 +178,11 @@ module material
homogenization_active homogenization_active
! BEGIN DEPRECATED ! BEGIN DEPRECATED
integer(pInt), dimension(:,:,:), allocatable, public :: phaseAt !< phase ID of every material point (ipc,ip,el) integer, dimension(:,:,:), allocatable, public :: phaseAt !< phase ID of every material point (ipc,ip,el)
integer(pInt), dimension(:,:,:), allocatable, public :: phasememberAt !< memberID of given phase at every material point (ipc,ip,el) integer, dimension(:,:,:), allocatable, public :: phasememberAt !< memberID of given phase at every material point (ipc,ip,el)
integer(pInt), dimension(:,:,:), allocatable, public, target :: mappingHomogenization !< mapping from material points to offset in heterogenous state/field integer, dimension(:,:,:), allocatable, public, target :: mappingHomogenization !< mapping from material points to offset in heterogenous state/field
integer(pInt), dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
! END DEPRECATED ! END DEPRECATED
type(tHomogMapping), allocatable, dimension(:), public :: & type(tHomogMapping), allocatable, dimension(:), public :: &
@ -256,13 +256,13 @@ subroutine material_init
use mesh, only: & use mesh, only: &
theMesh theMesh
integer(pInt), parameter :: FILEUNIT = 210_pInt integer, parameter :: FILEUNIT = 210
integer(pInt) :: m,c,h, myDebug, myPhase, myHomog integer :: m,c,h, myDebug, myPhase, myHomog
integer(pInt) :: & integer :: &
g, & !< grain number g, & !< grain number
i, & !< integration point number i, & !< integration point number
e !< element number e !< element number
integer(pInt), dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
CounterPhase, & CounterPhase, &
CounterHomogenization CounterHomogenization
@ -271,19 +271,19 @@ subroutine material_init
write(6,'(/,a)') ' <<<+- material init -+>>>' write(6,'(/,a)') ' <<<+- material init -+>>>'
call material_parsePhase() call material_parsePhase()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Phase parsed'; flush(6)
call material_parseMicrostructure() call material_parseMicrostructure()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Microstructure parsed'; flush(6)
call material_parseCrystallite() call material_parseCrystallite()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Crystallite parsed'; flush(6)
call material_parseHomogenization() call material_parseHomogenization()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Homogenization parsed'; flush(6)
call material_parseTexture() call material_parseTexture()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Texture parsed'; flush(6)
allocate(plasticState (size(config_phase))) allocate(plasticState (size(config_phase)))
allocate(sourceState (size(config_phase))) allocate(sourceState (size(config_phase)))
@ -303,34 +303,34 @@ subroutine material_init
allocate(temperatureRate (size(config_homogenization))) allocate(temperatureRate (size(config_homogenization)))
do m = 1_pInt,size(config_microstructure) do m = 1,size(config_microstructure)
if(microstructure_crystallite(m) < 1_pInt .or. & if(microstructure_crystallite(m) < 1 .or. &
microstructure_crystallite(m) > size(config_crystallite)) & microstructure_crystallite(m) > size(config_crystallite)) &
call IO_error(150_pInt,m,ext_msg='crystallite') call IO_error(150,m,ext_msg='crystallite')
if(minval(microstructure_phase(1:microstructure_Nconstituents(m),m)) < 1_pInt .or. & if(minval(microstructure_phase(1:microstructure_Nconstituents(m),m)) < 1 .or. &
maxval(microstructure_phase(1:microstructure_Nconstituents(m),m)) > size(config_phase)) & maxval(microstructure_phase(1:microstructure_Nconstituents(m),m)) > size(config_phase)) &
call IO_error(150_pInt,m,ext_msg='phase') call IO_error(150,m,ext_msg='phase')
if(minval(microstructure_texture(1:microstructure_Nconstituents(m),m)) < 1_pInt .or. & if(minval(microstructure_texture(1:microstructure_Nconstituents(m),m)) < 1 .or. &
maxval(microstructure_texture(1:microstructure_Nconstituents(m),m)) > size(config_texture)) & maxval(microstructure_texture(1:microstructure_Nconstituents(m),m)) > size(config_texture)) &
call IO_error(150_pInt,m,ext_msg='texture') call IO_error(150,m,ext_msg='texture')
if(microstructure_Nconstituents(m) < 1_pInt) & if(microstructure_Nconstituents(m) < 1) &
call IO_error(151_pInt,m) call IO_error(151,m)
enddo enddo
debugOut: if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then debugOut: if (iand(myDebug,debug_levelExtensive) /= 0) then
write(6,'(/,a,/)') ' MATERIAL configuration' write(6,'(/,a,/)') ' MATERIAL configuration'
write(6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains' write(6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains'
do h = 1_pInt,size(config_homogenization) do h = 1,size(config_homogenization)
write(6,'(1x,a32,1x,a16,1x,i6)') homogenization_name(h),homogenization_type(h),homogenization_Ngrains(h) write(6,'(1x,a32,1x,a16,1x,i6)') homogenization_name(h),homogenization_type(h),homogenization_Ngrains(h)
enddo enddo
write(6,'(/,a14,18x,1x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents','homogeneous' write(6,'(/,a14,18x,1x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents','homogeneous'
do m = 1_pInt,size(config_microstructure) do m = 1,size(config_microstructure)
write(6,'(1x,a32,1x,i11,1x,i12,1x,l13)') microstructure_name(m), & write(6,'(1x,a32,1x,i11,1x,i12,1x,l13)') microstructure_name(m), &
microstructure_crystallite(m), & microstructure_crystallite(m), &
microstructure_Nconstituents(m), & microstructure_Nconstituents(m), &
microstructure_elemhomo(m) microstructure_elemhomo(m)
if (microstructure_Nconstituents(m) > 0_pInt) then if (microstructure_Nconstituents(m) > 0) then
do c = 1_pInt,microstructure_Nconstituents(m) do c = 1,microstructure_Nconstituents(m)
write(6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(c,m)),& write(6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(c,m)),&
texture_name(microstructure_texture(c,m)),& texture_name(microstructure_texture(c,m)),&
microstructure_fraction(c,m) microstructure_fraction(c,m)
@ -383,23 +383,23 @@ subroutine material_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN DEPRECATED ! BEGIN DEPRECATED
allocate(phaseAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt) allocate(phaseAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0)
allocate(phasememberAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt) allocate(phasememberAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0)
allocate(mappingHomogenization (2, theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt) allocate(mappingHomogenization (2, theMesh%elem%nIPs,theMesh%Nelems),source=0)
allocate(mappingHomogenizationConst( theMesh%elem%nIPs,theMesh%Nelems),source=1_pInt) allocate(mappingHomogenizationConst( theMesh%elem%nIPs,theMesh%Nelems),source=1)
CounterHomogenization=0 CounterHomogenization=0
CounterPhase =0 CounterPhase =0
do e = 1_pInt,theMesh%Nelems do e = 1,theMesh%Nelems
myHomog = theMesh%homogenizationAt(e) myHomog = theMesh%homogenizationAt(e)
do i = 1_pInt, theMesh%elem%nIPs do i = 1, theMesh%elem%nIPs
CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1_pInt CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1
mappingHomogenization(1:2,i,e) = [CounterHomogenization(myHomog),huge(1)] mappingHomogenization(1:2,i,e) = [CounterHomogenization(myHomog),huge(1)]
do g = 1_pInt,homogenization_Ngrains(myHomog) do g = 1,homogenization_Ngrains(myHomog)
myPhase = material_phase(g,i,e) myPhase = material_phase(g,i,e)
CounterPhase(myPhase) = CounterPhase(myPhase)+1_pInt ! not distinguishing between instances of same phase CounterPhase(myPhase) = CounterPhase(myPhase)+1 ! not distinguishing between instances of same phase
phaseAt(g,i,e) = myPhase phaseAt(g,i,e) = myPhase
phasememberAt(g,i,e) = CounterPhase(myPhase) phasememberAt(g,i,e) = CounterPhase(myPhase)
enddo enddo
@ -429,33 +429,33 @@ subroutine material_parseHomogenization
use IO, only: & use IO, only: &
IO_error IO_error
integer(pInt) :: h integer :: h
character(len=65536) :: tag character(len=65536) :: tag
allocate(homogenization_type(size(config_homogenization)), source=HOMOGENIZATION_undefined_ID) allocate(homogenization_type(size(config_homogenization)), source=HOMOGENIZATION_undefined_ID)
allocate(thermal_type(size(config_homogenization)), source=THERMAL_isothermal_ID) allocate(thermal_type(size(config_homogenization)), source=THERMAL_isothermal_ID)
allocate(damage_type (size(config_homogenization)), source=DAMAGE_none_ID) allocate(damage_type (size(config_homogenization)), source=DAMAGE_none_ID)
allocate(homogenization_typeInstance(size(config_homogenization)), source=0_pInt) allocate(homogenization_typeInstance(size(config_homogenization)), source=0)
allocate(thermal_typeInstance(size(config_homogenization)), source=0_pInt) allocate(thermal_typeInstance(size(config_homogenization)), source=0)
allocate(damage_typeInstance(size(config_homogenization)), source=0_pInt) allocate(damage_typeInstance(size(config_homogenization)), source=0)
allocate(homogenization_Ngrains(size(config_homogenization)), source=0_pInt) allocate(homogenization_Ngrains(size(config_homogenization)), source=0)
allocate(homogenization_Noutput(size(config_homogenization)), source=0_pInt) allocate(homogenization_Noutput(size(config_homogenization)), source=0)
allocate(homogenization_active(size(config_homogenization)), source=.false.) !!!!!!!!!!!!!!! allocate(homogenization_active(size(config_homogenization)), source=.false.) !!!!!!!!!!!!!!!
allocate(thermal_initialT(size(config_homogenization)), source=300.0_pReal) allocate(thermal_initialT(size(config_homogenization)), source=300.0_pReal)
allocate(damage_initialPhi(size(config_homogenization)), source=1.0_pReal) allocate(damage_initialPhi(size(config_homogenization)), source=1.0_pReal)
forall (h = 1_pInt:size(config_homogenization)) & forall (h = 1:size(config_homogenization)) &
homogenization_active(h) = any(theMesh%homogenizationAt == h) homogenization_active(h) = any(theMesh%homogenizationAt == h)
do h=1_pInt, size(config_homogenization) do h=1, size(config_homogenization)
homogenization_Noutput(h) = config_homogenization(h)%countKeys('(output)') homogenization_Noutput(h) = config_homogenization(h)%countKeys('(output)')
tag = config_homogenization(h)%getString('mech') tag = config_homogenization(h)%getString('mech')
select case (trim(tag)) select case (trim(tag))
case(HOMOGENIZATION_NONE_label) case(HOMOGENIZATION_NONE_label)
homogenization_type(h) = HOMOGENIZATION_NONE_ID homogenization_type(h) = HOMOGENIZATION_NONE_ID
homogenization_Ngrains(h) = 1_pInt homogenization_Ngrains(h) = 1
case(HOMOGENIZATION_ISOSTRAIN_label) case(HOMOGENIZATION_ISOSTRAIN_label)
homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
homogenization_Ngrains(h) = config_homogenization(h)%getInt('nconstituents') homogenization_Ngrains(h) = config_homogenization(h)%getInt('nconstituents')
@ -463,7 +463,7 @@ subroutine material_parseHomogenization
homogenization_type(h) = HOMOGENIZATION_RGC_ID homogenization_type(h) = HOMOGENIZATION_RGC_ID
homogenization_Ngrains(h) = config_homogenization(h)%getInt('nconstituents') homogenization_Ngrains(h) = config_homogenization(h)%getInt('nconstituents')
case default case default
call IO_error(500_pInt,ext_msg=trim(tag)) call IO_error(500,ext_msg=trim(tag))
end select end select
homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h)) homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h))
@ -480,7 +480,7 @@ subroutine material_parseHomogenization
case(THERMAL_conduction_label) case(THERMAL_conduction_label)
thermal_type(h) = THERMAL_conduction_ID thermal_type(h) = THERMAL_conduction_ID
case default case default
call IO_error(500_pInt,ext_msg=trim(tag)) call IO_error(500,ext_msg=trim(tag))
end select end select
endif endif
@ -497,14 +497,14 @@ subroutine material_parseHomogenization
case(DAMAGE_NONLOCAL_label) case(DAMAGE_NONLOCAL_label)
damage_type(h) = DAMAGE_nonlocal_ID damage_type(h) = DAMAGE_nonlocal_ID
case default case default
call IO_error(500_pInt,ext_msg=trim(tag)) call IO_error(500,ext_msg=trim(tag))
end select end select
endif endif
enddo enddo
do h=1_pInt, size(config_homogenization) do h=1, size(config_homogenization)
homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h)) homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h))
thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h)) thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h))
damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h)) damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h))
@ -530,58 +530,58 @@ subroutine material_parseMicrostructure
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
strings strings
integer(pInt), allocatable, dimension(:) :: chunkPos integer, allocatable, dimension(:) :: chunkPos
integer(pInt) :: e, m, c, i integer :: e, m, c, i
character(len=65536) :: & character(len=65536) :: &
tag tag
allocate(microstructure_crystallite(size(config_microstructure)), source=0_pInt) allocate(microstructure_crystallite(size(config_microstructure)), source=0)
allocate(microstructure_Nconstituents(size(config_microstructure)), source=0_pInt) allocate(microstructure_Nconstituents(size(config_microstructure)), source=0)
allocate(microstructure_active(size(config_microstructure)), source=.false.) allocate(microstructure_active(size(config_microstructure)), source=.false.)
allocate(microstructure_elemhomo(size(config_microstructure)), source=.false.) allocate(microstructure_elemhomo(size(config_microstructure)), source=.false.)
if(any(theMesh%microstructureAt > size(config_microstructure))) & if(any(theMesh%microstructureAt > size(config_microstructure))) &
call IO_error(155_pInt,ext_msg='More microstructures in geometry than sections in material.config') call IO_error(155,ext_msg='More microstructures in geometry than sections in material.config')
forall (e = 1_pInt:theMesh%Nelems) & forall (e = 1:theMesh%Nelems) &
microstructure_active(theMesh%microstructureAt(e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements microstructure_active(theMesh%microstructureAt(e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements
do m=1_pInt, size(config_microstructure) do m=1, size(config_microstructure)
microstructure_Nconstituents(m) = config_microstructure(m)%countKeys('(constituent)') microstructure_Nconstituents(m) = config_microstructure(m)%countKeys('(constituent)')
microstructure_crystallite(m) = config_microstructure(m)%getInt('crystallite') microstructure_crystallite(m) = config_microstructure(m)%getInt('crystallite')
microstructure_elemhomo(m) = config_microstructure(m)%keyExists('/elementhomogeneous/') microstructure_elemhomo(m) = config_microstructure(m)%keyExists('/elementhomogeneous/')
enddo enddo
microstructure_maxNconstituents = maxval(microstructure_Nconstituents) microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
allocate(microstructure_phase (microstructure_maxNconstituents,size(config_microstructure)),source=0_pInt) allocate(microstructure_phase (microstructure_maxNconstituents,size(config_microstructure)),source=0)
allocate(microstructure_texture (microstructure_maxNconstituents,size(config_microstructure)),source=0_pInt) allocate(microstructure_texture (microstructure_maxNconstituents,size(config_microstructure)),source=0)
allocate(microstructure_fraction(microstructure_maxNconstituents,size(config_microstructure)),source=0.0_pReal) allocate(microstructure_fraction(microstructure_maxNconstituents,size(config_microstructure)),source=0.0_pReal)
allocate(strings(1)) ! Intel 16.0 Bug allocate(strings(1)) ! Intel 16.0 Bug
do m=1_pInt, size(config_microstructure) do m=1, size(config_microstructure)
strings = config_microstructure(m)%getStrings('(constituent)',raw=.true.) strings = config_microstructure(m)%getStrings('(constituent)',raw=.true.)
do c = 1_pInt, size(strings) do c = 1, size(strings)
chunkPos = IO_stringPos(strings(c)) chunkPos = IO_stringPos(strings(c))
do i = 1_pInt,5_pInt,2_pInt do i = 1,5,2
tag = IO_stringValue(strings(c),chunkPos,i) tag = IO_stringValue(strings(c),chunkPos,i)
select case (tag) select case (tag)
case('phase') case('phase')
microstructure_phase(c,m) = IO_intValue(strings(c),chunkPos,i+1_pInt) microstructure_phase(c,m) = IO_intValue(strings(c),chunkPos,i+1)
case('texture') case('texture')
microstructure_texture(c,m) = IO_intValue(strings(c),chunkPos,i+1_pInt) microstructure_texture(c,m) = IO_intValue(strings(c),chunkPos,i+1)
case('fraction') case('fraction')
microstructure_fraction(c,m) = IO_floatValue(strings(c),chunkPos,i+1_pInt) microstructure_fraction(c,m) = IO_floatValue(strings(c),chunkPos,i+1)
end select end select
enddo enddo
enddo enddo
enddo enddo
do m = 1_pInt, size(config_microstructure) do m = 1, size(config_microstructure)
if (dNeq(sum(microstructure_fraction(:,m)),1.0_pReal)) & if (dNeq(sum(microstructure_fraction(:,m)),1.0_pReal)) &
call IO_error(153_pInt,ext_msg=microstructure_name(m)) call IO_error(153,ext_msg=microstructure_name(m))
enddo enddo
end subroutine material_parseMicrostructure end subroutine material_parseMicrostructure
@ -592,10 +592,10 @@ end subroutine material_parseMicrostructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine material_parseCrystallite subroutine material_parseCrystallite
integer(pInt) :: c integer :: c
allocate(crystallite_Noutput(size(config_crystallite)),source=0_pInt) allocate(crystallite_Noutput(size(config_crystallite)),source=0)
do c=1_pInt, size(config_crystallite) do c=1, size(config_crystallite)
crystallite_Noutput(c) = config_crystallite(c)%countKeys('(output)') crystallite_Noutput(c) = config_crystallite(c)%countKeys('(output)')
enddo enddo
@ -611,19 +611,19 @@ subroutine material_parsePhase
IO_getTag, & IO_getTag, &
IO_stringValue IO_stringValue
integer(pInt) :: sourceCtr, kinematicsCtr, stiffDegradationCtr, p integer :: sourceCtr, kinematicsCtr, stiffDegradationCtr, p
character(len=65536), dimension(:), allocatable :: str character(len=65536), dimension(:), allocatable :: str
allocate(phase_elasticity(size(config_phase)),source=ELASTICITY_undefined_ID) allocate(phase_elasticity(size(config_phase)),source=ELASTICITY_undefined_ID)
allocate(phase_plasticity(size(config_phase)),source=PLASTICITY_undefined_ID) allocate(phase_plasticity(size(config_phase)),source=PLASTICITY_undefined_ID)
allocate(phase_Nsources(size(config_phase)), source=0_pInt) allocate(phase_Nsources(size(config_phase)), source=0)
allocate(phase_Nkinematics(size(config_phase)), source=0_pInt) allocate(phase_Nkinematics(size(config_phase)), source=0)
allocate(phase_NstiffnessDegradations(size(config_phase)),source=0_pInt) allocate(phase_NstiffnessDegradations(size(config_phase)),source=0)
allocate(phase_Noutput(size(config_phase)), source=0_pInt) allocate(phase_Noutput(size(config_phase)), source=0)
allocate(phase_localPlasticity(size(config_phase)), source=.false.) allocate(phase_localPlasticity(size(config_phase)), source=.false.)
do p=1_pInt, size(config_phase) do p=1, size(config_phase)
phase_Noutput(p) = config_phase(p)%countKeys('(output)') phase_Noutput(p) = config_phase(p)%countKeys('(output)')
phase_Nsources(p) = config_phase(p)%countKeys('(source)') phase_Nsources(p) = config_phase(p)%countKeys('(source)')
phase_Nkinematics(p) = config_phase(p)%countKeys('(kinematics)') phase_Nkinematics(p) = config_phase(p)%countKeys('(kinematics)')
@ -634,7 +634,7 @@ subroutine material_parsePhase
case (ELASTICITY_HOOKE_label) case (ELASTICITY_HOOKE_label)
phase_elasticity(p) = ELASTICITY_HOOKE_ID phase_elasticity(p) = ELASTICITY_HOOKE_ID
case default case default
call IO_error(200_pInt,ext_msg=trim(config_phase(p)%getString('elasticity'))) call IO_error(200,ext_msg=trim(config_phase(p)%getString('elasticity')))
end select end select
select case (config_phase(p)%getString('plasticity')) select case (config_phase(p)%getString('plasticity'))
@ -653,7 +653,7 @@ subroutine material_parsePhase
case (PLASTICITY_NONLOCAL_label) case (PLASTICITY_NONLOCAL_label)
phase_plasticity(p) = PLASTICITY_NONLOCAL_ID phase_plasticity(p) = PLASTICITY_NONLOCAL_ID
case default case default
call IO_error(201_pInt,ext_msg=trim(config_phase(p)%getString('plasticity'))) call IO_error(201,ext_msg=trim(config_phase(p)%getString('plasticity')))
end select end select
enddo enddo
@ -662,7 +662,7 @@ subroutine material_parsePhase
allocate(phase_kinematics(maxval(phase_Nkinematics),size(config_phase)), source=KINEMATICS_undefined_ID) allocate(phase_kinematics(maxval(phase_Nkinematics),size(config_phase)), source=KINEMATICS_undefined_ID)
allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),size(config_phase)), & allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),size(config_phase)), &
source=STIFFNESS_DEGRADATION_undefined_ID) source=STIFFNESS_DEGRADATION_undefined_ID)
do p=1_pInt, size(config_phase) do p=1, size(config_phase)
#if defined(__GFORTRAN__) || defined(__PGI) #if defined(__GFORTRAN__) || defined(__PGI)
str = ['GfortranBug86277'] str = ['GfortranBug86277']
str = config_phase(p)%getStrings('(source)',defaultVal=str) str = config_phase(p)%getStrings('(source)',defaultVal=str)
@ -670,7 +670,7 @@ subroutine material_parsePhase
#else #else
str = config_phase(p)%getStrings('(source)',defaultVal=[character(len=65536)::]) str = config_phase(p)%getStrings('(source)',defaultVal=[character(len=65536)::])
#endif #endif
do sourceCtr = 1_pInt, size(str) do sourceCtr = 1, size(str)
select case (trim(str(sourceCtr))) select case (trim(str(sourceCtr)))
case (SOURCE_thermal_dissipation_label) case (SOURCE_thermal_dissipation_label)
phase_source(sourceCtr,p) = SOURCE_thermal_dissipation_ID phase_source(sourceCtr,p) = SOURCE_thermal_dissipation_ID
@ -694,7 +694,7 @@ subroutine material_parsePhase
#else #else
str = config_phase(p)%getStrings('(kinematics)',defaultVal=[character(len=65536)::]) str = config_phase(p)%getStrings('(kinematics)',defaultVal=[character(len=65536)::])
#endif #endif
do kinematicsCtr = 1_pInt, size(str) do kinematicsCtr = 1, size(str)
select case (trim(str(kinematicsCtr))) select case (trim(str(kinematicsCtr)))
case (KINEMATICS_cleavage_opening_label) case (KINEMATICS_cleavage_opening_label)
phase_kinematics(kinematicsCtr,p) = KINEMATICS_cleavage_opening_ID phase_kinematics(kinematicsCtr,p) = KINEMATICS_cleavage_opening_ID
@ -711,7 +711,7 @@ subroutine material_parsePhase
#else #else
str = config_phase(p)%getStrings('(stiffness_degradation)',defaultVal=[character(len=65536)::]) str = config_phase(p)%getStrings('(stiffness_degradation)',defaultVal=[character(len=65536)::])
#endif #endif
do stiffDegradationCtr = 1_pInt, size(str) do stiffDegradationCtr = 1, size(str)
select case (trim(str(stiffDegradationCtr))) select case (trim(str(stiffDegradationCtr)))
case (STIFFNESS_DEGRADATION_damage_label) case (STIFFNESS_DEGRADATION_damage_label)
phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID
@ -719,10 +719,10 @@ subroutine material_parsePhase
enddo enddo
enddo enddo
allocate(phase_plasticityInstance(size(config_phase)), source=0_pInt) allocate(phase_plasticityInstance(size(config_phase)), source=0)
allocate(phase_elasticityInstance(size(config_phase)), source=0_pInt) allocate(phase_elasticityInstance(size(config_phase)), source=0)
do p=1_pInt, size(config_phase) do p=1, size(config_phase)
phase_elasticityInstance(p) = count(phase_elasticity(1:p) == phase_elasticity(p)) phase_elasticityInstance(p) = count(phase_elasticity(1:p) == phase_elasticity(p))
phase_plasticityInstance(p) = count(phase_plasticity(1:p) == phase_plasticity(p)) phase_plasticityInstance(p) = count(phase_plasticity(1:p) == phase_plasticity(p))
enddo enddo
@ -739,13 +739,13 @@ subroutine material_parseTexture
IO_floatValue, & IO_floatValue, &
IO_stringValue IO_stringValue
integer(pInt) :: section, gauss, j, t, i integer :: section, gauss, j, t, i
character(len=65536), dimension(:), allocatable :: strings ! Values for given key in material config character(len=65536), dimension(:), allocatable :: strings ! Values for given key in material config
integer(pInt), dimension(:), allocatable :: chunkPos integer, dimension(:), allocatable :: chunkPos
allocate(texture_Ngauss(size(config_texture)), source=0_pInt) allocate(texture_Ngauss(size(config_texture)), source=0)
do t=1_pInt, size(config_texture) do t=1, size(config_texture)
texture_Ngauss(t) = config_texture(t)%countKeys('(gauss)') texture_Ngauss(t) = config_texture(t)%countKeys('(gauss)')
if (config_texture(t)%keyExists('symmetry')) call IO_error(147,ext_msg='symmetry') if (config_texture(t)%keyExists('symmetry')) call IO_error(147,ext_msg='symmetry')
if (config_texture(t)%keyExists('(random)')) call IO_error(147,ext_msg='(random)') if (config_texture(t)%keyExists('(random)')) call IO_error(147,ext_msg='(random)')
@ -757,13 +757,13 @@ subroutine material_parseTexture
allocate(texture_transformation(3,3,size(config_texture)), source=0.0_pReal) allocate(texture_transformation(3,3,size(config_texture)), source=0.0_pReal)
texture_transformation = spread(math_I3,3,size(config_texture)) texture_transformation = spread(math_I3,3,size(config_texture))
do t=1_pInt, size(config_texture) do t=1, size(config_texture)
section = t section = t
gauss = 0_pInt gauss = 0
if (config_texture(t)%keyExists('axes')) then if (config_texture(t)%keyExists('axes')) then
strings = config_texture(t)%getStrings('axes') strings = config_texture(t)%getStrings('axes')
do j = 1_pInt, 3_pInt ! look for "x", "y", and "z" entries do j = 1, 3 ! look for "x", "y", and "z" entries
select case (strings(j)) select case (strings(j))
case('x', '+x') case('x', '+x')
texture_transformation(j,1:3,t) = [ 1.0_pReal, 0.0_pReal, 0.0_pReal] ! original axis is now +x-axis texture_transformation(j,1:3,t) = [ 1.0_pReal, 0.0_pReal, 0.0_pReal] ! original axis is now +x-axis
@ -778,25 +778,25 @@ subroutine material_parseTexture
case('-z') case('-z')
texture_transformation(j,1:3,t) = [ 0.0_pReal, 0.0_pReal,-1.0_pReal] ! original axis is now -z-axis texture_transformation(j,1:3,t) = [ 0.0_pReal, 0.0_pReal,-1.0_pReal] ! original axis is now -z-axis
case default case default
call IO_error(157_pInt,t) call IO_error(157,t)
end select end select
enddo enddo
if(dNeq(math_det33(texture_transformation(1:3,1:3,t)),1.0_pReal)) call IO_error(157_pInt,t) if(dNeq(math_det33(texture_transformation(1:3,1:3,t)),1.0_pReal)) call IO_error(157,t)
endif endif
if (config_texture(t)%keyExists('(gauss)')) then if (config_texture(t)%keyExists('(gauss)')) then
gauss = gauss + 1_pInt gauss = gauss + 1
strings = config_texture(t)%getStrings('(gauss)',raw= .true.) strings = config_texture(t)%getStrings('(gauss)',raw= .true.)
do i = 1_pInt , size(strings) do i = 1 , size(strings)
chunkPos = IO_stringPos(strings(i)) chunkPos = IO_stringPos(strings(i))
do j = 1_pInt,9_pInt,2_pInt do j = 1,9,2
select case (IO_stringValue(strings(i),chunkPos,j)) select case (IO_stringValue(strings(i),chunkPos,j))
case('phi1') case('phi1')
texture_Gauss(1,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad texture_Gauss(1,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
case('phi') case('phi')
texture_Gauss(2,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad texture_Gauss(2,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
case('phi2') case('phi2')
texture_Gauss(3,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad texture_Gauss(3,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
end select end select
enddo enddo
enddo enddo
@ -817,7 +817,7 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,&
use numerics, only: & use numerics, only: &
numerics_integrator numerics_integrator
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
NofMyPhase, & NofMyPhase, &
sizeState, & sizeState, &
@ -842,13 +842,13 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,&
allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 1_pInt) then if (numerics_integrator == 1) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif endif
if (numerics_integrator == 4_pInt) & if (numerics_integrator == 4) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 5_pInt) & if (numerics_integrator == 5) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
@ -864,7 +864,7 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,&
use numerics, only: & use numerics, only: &
numerics_integrator numerics_integrator
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
of, & of, &
NofMyPhase, & NofMyPhase, &
@ -882,13 +882,13 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,&
allocate(sourceState(phase)%p(of)%state (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 1_pInt) then if (numerics_integrator == 1) then
allocate(sourceState(phase)%p(of)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif endif
if (numerics_integrator == 4_pInt) & if (numerics_integrator == 4) &
allocate(sourceState(phase)%p(of)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 5_pInt) & if (numerics_integrator == 5) &
allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
@ -905,10 +905,10 @@ subroutine material_populateGrains
use mesh, only: & use mesh, only: &
theMesh theMesh
integer(pInt) :: e,i,c,homog,micro integer :: e,i,c,homog,micro
allocate(material_phase(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0_pInt) allocate(material_phase(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0)
allocate(material_texture(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0_pInt) allocate(material_texture(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0)
allocate(material_EulerAngles(3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0.0_pReal) allocate(material_EulerAngles(3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0.0_pReal)
do e = 1, theMesh%Nelems do e = 1, theMesh%Nelems

971
src/mesh_abaqus.f90
View File

File diff suppressed because it is too large Load Diff

750
src/mesh_marc.f90
View File

File diff suppressed because it is too large Load Diff

78
src/source_damage_anisoBrittle.f90
View File

@ -5,23 +5,21 @@
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_damage_anisoBrittle module source_damage_anisoBrittle
use prec, only: & use prec
pReal, &
pInt
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
source_damage_anisoBrittle_instance !< instance of source mechanism source_damage_anisoBrittle_instance !< instance of source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: & integer, dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_sizePostResult !< size of each post result output source_damage_anisoBrittle_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_output !< name of each post result output source_damage_anisoBrittle_output !< name of each post result output
integer(pInt), dimension(:,:), allocatable, private :: & integer, dimension(:,:), allocatable, private :: &
source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family
enum, bind(c) enum, bind(c)
@ -40,9 +38,9 @@ module source_damage_anisoBrittle
critLoad critLoad
real(pReal), dimension(:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems cleavage_systems
integer(pInt) :: & integer :: &
totalNcleavage totalNcleavage
integer(pInt), dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
Ncleavage Ncleavage
integer(kind(undefined_ID)), allocatable, dimension(:) :: & integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output outputID !< ID of each post result output
@ -65,8 +63,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_init subroutine source_damage_anisoBrittle_init
use prec, only: &
pStringLen
use debug, only: & use debug, only: &
debug_level,& debug_level,&
debug_constitutive,& debug_constitutive,&
@ -91,10 +87,10 @@ subroutine source_damage_anisoBrittle_init
lattice_SchmidMatrix_cleavage, & lattice_SchmidMatrix_cleavage, &
lattice_maxNcleavageFamily lattice_maxNcleavageFamily
integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer :: Ninstance,phase,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p ,i integer :: NofMyPhase,p ,i
integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] integer, dimension(0), parameter :: emptyIntArray = [integer::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: & integer(kind(undefined_ID)) :: &
outputID outputID
@ -105,14 +101,14 @@ subroutine source_damage_anisoBrittle_init
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>' write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'
Ninstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt) Ninstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID))
if (Ninstance == 0_pInt) return if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0_pInt) allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0)
allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0_pInt) allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0)
do phase = 1, material_Nphase do phase = 1, material_Nphase
source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID) source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)
@ -121,11 +117,11 @@ subroutine source_damage_anisoBrittle_init
enddo enddo
enddo enddo
allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0_pInt) allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0)
allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance)) allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoBrittle_output = '' source_damage_anisoBrittle_output = ''
allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0_pInt) allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0)
allocate(param(Ninstance)) allocate(param(Ninstance))
@ -162,18 +158,18 @@ subroutine source_damage_anisoBrittle_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range ! exit if any parameter is out of range
if (extmsg /= '') & if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output pararameters ! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray) outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0)) allocate(prm%outputID(0))
do i=1_pInt, size(outputs) do i=1, size(outputs)
outputID = undefined_ID outputID = undefined_ID
select case(outputs(i)) select case(outputs(i))
case ('anisobrittle_drivingforce') case ('anisobrittle_drivingforce')
source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1_pInt source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1
source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i) source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID] prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -189,7 +185,7 @@ subroutine source_damage_anisoBrittle_init
sourceOffset = source_damage_anisoBrittle_offset(phase) sourceOffset = source_damage_anisoBrittle_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt) call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance)) sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
@ -217,13 +213,13 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
lattice_maxNcleavageFamily, & lattice_maxNcleavageFamily, &
lattice_NcleavageSystem lattice_NcleavageSystem
integer(pInt), intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S S
integer(pInt) :: & integer :: &
phase, & phase, &
constituent, & constituent, &
instance, & instance, &
@ -243,10 +239,10 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
index = 1_pInt index = 1
do f = 1_pInt,lattice_maxNcleavageFamily do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1_pInt,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family do i = 1,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
@ -263,7 +259,7 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ & (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ &
param(instance)%critDisp(index) param(instance)%critDisp(index)
index = index + 1_pInt index = index + 1
enddo enddo
enddo enddo
@ -276,7 +272,7 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
use material, only: & use material, only: &
sourceState sourceState
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
constituent constituent
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -284,7 +280,7 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer(pInt) :: & integer :: &
sourceOffset sourceOffset
sourceOffset = source_damage_anisoBrittle_offset(phase) sourceOffset = source_damage_anisoBrittle_offset(phase)
@ -303,27 +299,27 @@ function source_damage_anisoBrittle_postResults(phase, constituent)
use material, only: & use material, only: &
sourceState sourceState
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
constituent constituent
real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, & real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, &
source_damage_anisoBrittle_instance(phase)))) :: & source_damage_anisoBrittle_instance(phase)))) :: &
source_damage_anisoBrittle_postResults source_damage_anisoBrittle_postResults
integer(pInt) :: & integer :: &
instance, sourceOffset, o, c instance, sourceOffset, o, c
instance = source_damage_anisoBrittle_instance(phase) instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase) sourceOffset = source_damage_anisoBrittle_offset(phase)
c = 0_pInt c = 0
do o = 1_pInt,size(param(instance)%outputID) do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o)) select case(param(instance)%outputID(o))
case (damage_drivingforce_ID) case (damage_drivingforce_ID)
source_damage_anisoBrittle_postResults(c+1_pInt) = & source_damage_anisoBrittle_postResults(c+1) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent) sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1_pInt c = c + 1
end select end select
enddo enddo

58
src/source_damage_anisoDuctile.f90
View File

@ -5,20 +5,18 @@
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_damage_anisoDuctile module source_damage_anisoDuctile
use prec, only: & use prec
pReal, &
pInt
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism source_damage_anisoDuctile_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: & integer, dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_sizePostResult !< size of each post result output source_damage_anisoDuctile_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_output !< name of each post result output source_damage_anisoDuctile_output !< name of each post result output
@ -59,8 +57,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_init subroutine source_damage_anisoDuctile_init
use prec, only: &
pStringLen
use debug, only: & use debug, only: &
debug_level,& debug_level,&
debug_constitutive,& debug_constitutive,&
@ -82,11 +78,11 @@ subroutine source_damage_anisoDuctile_init
config_phase config_phase
integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer :: Ninstance,phase,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p ,i integer :: NofMyPhase,p ,i
integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] integer, dimension(0), parameter :: emptyIntArray = [integer::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: & integer(kind(undefined_ID)) :: &
outputID outputID
@ -98,13 +94,13 @@ subroutine source_damage_anisoDuctile_init
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>' write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'
Ninstance = count(phase_source == SOURCE_damage_anisoDuctile_ID) Ninstance = count(phase_source == SOURCE_damage_anisoDuctile_ID)
if (Ninstance == 0_pInt) return if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoDuctile_offset(size(config_phase)), source=0_pInt) allocate(source_damage_anisoDuctile_offset(size(config_phase)), source=0)
allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0_pInt) allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0)
do phase = 1, size(config_phase) do phase = 1, size(config_phase)
source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID) source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)
@ -113,7 +109,7 @@ subroutine source_damage_anisoDuctile_init
enddo enddo
enddo enddo
allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance)) allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoDuctile_output = '' source_damage_anisoDuctile_output = ''
@ -146,18 +142,18 @@ subroutine source_damage_anisoDuctile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range ! exit if any parameter is out of range
if (extmsg /= '') & if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output pararameters ! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray) outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0)) allocate(prm%outputID(0))
do i=1_pInt, size(outputs) do i=1, size(outputs)
outputID = undefined_ID outputID = undefined_ID
select case(outputs(i)) select case(outputs(i))
case ('anisoductile_drivingforce') case ('anisoductile_drivingforce')
source_damage_anisoDuctile_sizePostResult(i,source_damage_anisoDuctile_instance(p)) = 1_pInt source_damage_anisoDuctile_sizePostResult(i,source_damage_anisoDuctile_instance(p)) = 1
source_damage_anisoDuctile_output(i,source_damage_anisoDuctile_instance(p)) = outputs(i) source_damage_anisoDuctile_output(i,source_damage_anisoDuctile_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID] prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -173,7 +169,7 @@ subroutine source_damage_anisoDuctile_init
instance = source_damage_anisoDuctile_instance(phase) instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase) sourceOffset = source_damage_anisoDuctile_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt) call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoDuctile_sizePostResult(:,instance)) sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoDuctile_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
@ -193,11 +189,11 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
damage, & damage, &
damageMapping damageMapping
integer(pInt), intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
integer(pInt) :: & integer :: &
phase, & phase, &
constituent, & constituent, &
sourceOffset, & sourceOffset, &
@ -229,7 +225,7 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph
use material, only: & use material, only: &
sourceState sourceState
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
constituent constituent
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -237,7 +233,7 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer(pInt) :: & integer :: &
sourceOffset sourceOffset
sourceOffset = source_damage_anisoDuctile_offset(phase) sourceOffset = source_damage_anisoDuctile_offset(phase)
@ -256,27 +252,27 @@ function source_damage_anisoDuctile_postResults(phase, constituent)
use material, only: & use material, only: &
sourceState sourceState
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
constituent constituent
real(pReal), dimension(sum(source_damage_anisoDuctile_sizePostResult(:, & real(pReal), dimension(sum(source_damage_anisoDuctile_sizePostResult(:, &
source_damage_anisoDuctile_instance(phase)))) :: & source_damage_anisoDuctile_instance(phase)))) :: &
source_damage_anisoDuctile_postResults source_damage_anisoDuctile_postResults
integer(pInt) :: & integer :: &
instance, sourceOffset, o, c instance, sourceOffset, o, c
instance = source_damage_anisoDuctile_instance(phase) instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase) sourceOffset = source_damage_anisoDuctile_offset(phase)
c = 0_pInt c = 0
do o = 1_pInt,size(param(instance)%outputID) do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o)) select case(param(instance)%outputID(o))
case (damage_drivingforce_ID) case (damage_drivingforce_ID)
source_damage_anisoDuctile_postResults(c+1_pInt) = & source_damage_anisoDuctile_postResults(c+1) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent) sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1_pInt c = c + 1
end select end select
enddo enddo

56
src/source_damage_isoBrittle.f90
View File

@ -5,20 +5,18 @@
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_damage_isoBrittle module source_damage_isoBrittle
use prec, only: & use prec
pReal, &
pInt
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
source_damage_isoBrittle_offset, & !< which source is my current damage mechanism? source_damage_isoBrittle_offset, & !< which source is my current damage mechanism?
source_damage_isoBrittle_instance !< instance of damage source mechanism source_damage_isoBrittle_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: & integer, dimension(:,:), allocatable, target, public :: &
source_damage_isoBrittle_sizePostResult !< size of each post result output source_damage_isoBrittle_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_isoBrittle_output !< name of each post result output source_damage_isoBrittle_output !< name of each post result output
enum, bind(c) enum, bind(c)
@ -53,8 +51,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_init subroutine source_damage_isoBrittle_init
use prec, only: &
pStringLen
use debug, only: & use debug, only: &
debug_level,& debug_level,&
debug_constitutive,& debug_constitutive,&
@ -75,9 +71,9 @@ subroutine source_damage_isoBrittle_init
material_Nphase material_Nphase
integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer :: Ninstance,phase,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p,i integer :: NofMyPhase,p,i
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: & integer(kind(undefined_ID)) :: &
outputID outputID
@ -88,14 +84,14 @@ subroutine source_damage_isoBrittle_init
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>' write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>'
Ninstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt) Ninstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID))
if (Ninstance == 0_pInt) return if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_isoBrittle_offset(material_Nphase), source=0_pInt) allocate(source_damage_isoBrittle_offset(material_Nphase), source=0)
allocate(source_damage_isoBrittle_instance(material_Nphase), source=0_pInt) allocate(source_damage_isoBrittle_instance(material_Nphase), source=0)
do phase = 1, material_Nphase do phase = 1, material_Nphase
source_damage_isoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoBrittle_ID) source_damage_isoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoBrittle_ID)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)
@ -104,7 +100,7 @@ subroutine source_damage_isoBrittle_init
enddo enddo
enddo enddo
allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),Ninstance)) allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),Ninstance))
source_damage_isoBrittle_output = '' source_damage_isoBrittle_output = ''
@ -129,18 +125,18 @@ subroutine source_damage_isoBrittle_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range ! exit if any parameter is out of range
if (extmsg /= '') & if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISOBRITTLE_LABEL//')') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISOBRITTLE_LABEL//')')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output pararameters ! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray) outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0)) allocate(prm%outputID(0))
do i=1_pInt, size(outputs) do i=1, size(outputs)
outputID = undefined_ID outputID = undefined_ID
select case(outputs(i)) select case(outputs(i))
case ('isobrittle_drivingforce') case ('isobrittle_drivingforce')
source_damage_isoBrittle_sizePostResult(i,source_damage_isoBrittle_instance(p)) = 1_pInt source_damage_isoBrittle_sizePostResult(i,source_damage_isoBrittle_instance(p)) = 1
source_damage_isoBrittle_output(i,source_damage_isoBrittle_instance(p)) = outputs(i) source_damage_isoBrittle_output(i,source_damage_isoBrittle_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID] prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -156,7 +152,7 @@ subroutine source_damage_isoBrittle_init
instance = source_damage_isoBrittle_instance(phase) instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase) sourceOffset = source_damage_isoBrittle_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,1_pInt) call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,1)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoBrittle_sizePostResult(:,instance)) sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoBrittle_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
@ -175,7 +171,7 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
math_sym33to6, & math_sym33to6, &
math_I3 math_I3
integer(pInt), intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
@ -183,7 +179,7 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
Fe Fe
real(pReal), intent(in), dimension(6,6) :: & real(pReal), intent(in), dimension(6,6) :: &
C C
integer(pInt) :: & integer :: &
phase, constituent, instance, sourceOffset phase, constituent, instance, sourceOffset
real(pReal) :: & real(pReal) :: &
strain(6), & strain(6), &
@ -219,7 +215,7 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD
use material, only: & use material, only: &
sourceState sourceState
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
constituent constituent
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -227,7 +223,7 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer(pInt) :: & integer :: &
instance, sourceOffset instance, sourceOffset
instance = source_damage_isoBrittle_instance(phase) instance = source_damage_isoBrittle_instance(phase)
@ -248,25 +244,25 @@ function source_damage_isoBrittle_postResults(phase, constituent)
use material, only: & use material, only: &
sourceState sourceState
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
constituent constituent
real(pReal), dimension(sum(source_damage_isoBrittle_sizePostResult(:, & real(pReal), dimension(sum(source_damage_isoBrittle_sizePostResult(:, &
source_damage_isoBrittle_instance(phase)))) :: & source_damage_isoBrittle_instance(phase)))) :: &
source_damage_isoBrittle_postResults source_damage_isoBrittle_postResults
integer(pInt) :: & integer :: &
instance, sourceOffset, o, c instance, sourceOffset, o, c
instance = source_damage_isoBrittle_instance(phase) instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase) sourceOffset = source_damage_isoBrittle_offset(phase)
c = 0_pInt c = 0
do o = 1_pInt,size(param(instance)%outputID) do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o)) select case(param(instance)%outputID(o))
case (damage_drivingforce_ID) case (damage_drivingforce_ID)
source_damage_isoBrittle_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent) source_damage_isoBrittle_postResults(c+1) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1 c = c + 1
end select end select

50
src/source_damage_isoDuctile.f90
View File

@ -5,20 +5,18 @@
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_damage_isoDuctile module source_damage_isoDuctile
use prec, only: & use prec
pReal, &
pInt
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_isoDuctile_instance !< instance of damage source mechanism source_damage_isoDuctile_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: & integer, dimension(:,:), allocatable, target, public :: &
source_damage_isoDuctile_sizePostResult !< size of each post result output source_damage_isoDuctile_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_isoDuctile_output !< name of each post result output source_damage_isoDuctile_output !< name of each post result output
@ -53,8 +51,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_init subroutine source_damage_isoDuctile_init
use prec, only: &
pStringLen
use debug, only: & use debug, only: &
debug_level,& debug_level,&
debug_constitutive,& debug_constitutive,&
@ -75,8 +71,8 @@ subroutine source_damage_isoDuctile_init
material_Nphase material_Nphase
integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer :: Ninstance,phase,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p,i integer :: NofMyPhase,p,i
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: & integer(kind(undefined_ID)) :: &
outputID outputID
@ -89,13 +85,13 @@ subroutine source_damage_isoDuctile_init
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init -+>>>' write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init -+>>>'
Ninstance = count(phase_source == SOURCE_damage_isoDuctile_ID) Ninstance = count(phase_source == SOURCE_damage_isoDuctile_ID)
if (Ninstance == 0_pInt) return if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_isoDuctile_offset(material_Nphase), source=0_pInt) allocate(source_damage_isoDuctile_offset(material_Nphase), source=0)
allocate(source_damage_isoDuctile_instance(material_Nphase), source=0_pInt) allocate(source_damage_isoDuctile_instance(material_Nphase), source=0)
do phase = 1, material_Nphase do phase = 1, material_Nphase
source_damage_isoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoDuctile_ID) source_damage_isoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoDuctile_ID)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)
@ -104,7 +100,7 @@ subroutine source_damage_isoDuctile_init
enddo enddo
enddo enddo
allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),Ninstance)) allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),Ninstance))
source_damage_isoDuctile_output = '' source_damage_isoDuctile_output = ''
@ -129,18 +125,18 @@ subroutine source_damage_isoDuctile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range ! exit if any parameter is out of range
if (extmsg /= '') & if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISODUCTILE_LABEL//')') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISODUCTILE_LABEL//')')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output pararameters ! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray) outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0)) allocate(prm%outputID(0))
do i=1_pInt, size(outputs) do i=1, size(outputs)
outputID = undefined_ID outputID = undefined_ID
select case(outputs(i)) select case(outputs(i))
case ('isoductile_drivingforce') case ('isoductile_drivingforce')
source_damage_isoDuctile_sizePostResult(i,source_damage_isoDuctile_instance(p)) = 1_pInt source_damage_isoDuctile_sizePostResult(i,source_damage_isoDuctile_instance(p)) = 1
source_damage_isoDuctile_output(i,source_damage_isoDuctile_instance(p)) = outputs(i) source_damage_isoDuctile_output(i,source_damage_isoDuctile_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID] prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -155,7 +151,7 @@ subroutine source_damage_isoDuctile_init
instance = source_damage_isoDuctile_instance(phase) instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase) sourceOffset = source_damage_isoDuctile_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt) call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoDuctile_sizePostResult(:,instance)) sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoDuctile_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
@ -176,11 +172,11 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
damage, & damage, &
damageMapping damageMapping
integer(pInt), intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
integer(pInt) :: & integer :: &
phase, constituent, instance, homog, sourceOffset, damageOffset phase, constituent, instance, homog, sourceOffset, damageOffset
phase = phaseAt(ipc,ip,el) phase = phaseAt(ipc,ip,el)
@ -204,7 +200,7 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD
use material, only: & use material, only: &
sourceState sourceState
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
constituent constituent
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -212,7 +208,7 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer(pInt) :: & integer :: &
sourceOffset sourceOffset
sourceOffset = source_damage_isoDuctile_offset(phase) sourceOffset = source_damage_isoDuctile_offset(phase)
@ -231,25 +227,25 @@ function source_damage_isoDuctile_postResults(phase, constituent)
use material, only: & use material, only: &
sourceState sourceState
integer(pInt), intent(in) :: & integer, intent(in) :: &
phase, & phase, &
constituent constituent
real(pReal), dimension(sum(source_damage_isoDuctile_sizePostResult(:, & real(pReal), dimension(sum(source_damage_isoDuctile_sizePostResult(:, &
source_damage_isoDuctile_instance(phase)))) :: & source_damage_isoDuctile_instance(phase)))) :: &
source_damage_isoDuctile_postResults source_damage_isoDuctile_postResults
integer(pInt) :: & integer :: &
instance, sourceOffset, o, c instance, sourceOffset, o, c
instance = source_damage_isoDuctile_instance(phase) instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase) sourceOffset = source_damage_isoDuctile_offset(phase)
c = 0_pInt c = 0
do o = 1_pInt,size(param(instance)%outputID) do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o)) select case(param(instance)%outputID(o))
case (damage_drivingforce_ID) case (damage_drivingforce_ID)
source_damage_isoDuctile_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent) source_damage_isoDuctile_postResults(c+1) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1 c = c + 1
end select end select