revert back to single level design for source code

This commit is contained in:
zhangc43 2016-02-29 14:20:52 -05:00
parent e3f82f7629
commit ad6432c173
28 changed files with 32 additions and 18860 deletions

View File

@ -79,7 +79,7 @@ endif (DAMASK_DRIVER STREQUAL "SPECTRAL")
# set system include directories
include_directories(${CMAKE_SOURCE_DIR}
include_directories(${CMAKE_SOURCE_DIR}/code
${PETSC_DIR}/include
lib
${HDF5_DIR}/include

View File

@ -46,14 +46,38 @@ target_link_libraries(DAMASK_LATTICE DAMASK_MATERIAL)
add_library(DAMASK_DRIVERS ALIAS DAMASK_LATTICE)
# For each modular section
add_subdirectory(plastic)
# add_subdirectory(kinematics)
# add_subdirectory(source)
add_library (DAMASK_PLASTIC "plastic_dislotwin.f90"
"plastic_disloUCLA.f90"
"plastic_isotropic.f90"
"plastic_j2.f90"
"plastic_phenopowerlaw.f90"
"plastic_titanmod.f90"
"plastic_nonlocal.f90"
"plastic_none.f90"
"plastic_phenoplus.f90")
target_link_libraries(DAMASK_PLASTIC DAMASK_DRIVERS)
# add_library(DAMASK_CONSTITUTIVE "constitutive.f90")
# target_link_libraries(DAMASK_CONSTITUTIVE DAMASK_DRIVERS)
# add_library(DAMASK_DRIVERS ALIAS DAMASK_CONSTITUTIVE)
add_library (DAMASK_KINEMATICS "kinematics_cleavage_opening.f90"
"kinematics_slipplane_opening.f90"
"kinematics_thermal_expansion.f90"
"kinematics_vacancy_strain.f90"
"kinematics_hydrogen_strain.f90")
target_link_libraries(DAMASK_KINEMATICS DAMASK_DRIVERS)
add_library (DAMASK_SOURCE "source_thermal_dissipation.f90"
"source_thermal_externalheat.f90"
"source_damage_isoBrittle.f90"
"source_damage_isoDuctile.f90"
"source_damage_anisoBrittle.f90"
"source_damage_anisoDuctile.f90"
"source_vacancy_phenoplasticity.f90"
"source_vacancy_irradiation.f90"
"source_vacancy_thermalfluc.f90")
target_link_libraries(DAMASK_SOURCE DAMASK_DRIVERS)
add_library(DAMASK_CONSTITUTIVE "constitutive.f90")
target_link_libraries(DAMASK_CONSTITUTIVE DAMASK_PLASTIC DAMASK_KINEMATICS DAMASK_SOURCE)
# compile spectral solver
add_executable(DAMASKSpectral.exe DAMASK_spectral.f90)
target_link_libraries (DAMASKSpectral.exe DAMASK_PLASTIC)
target_link_libraries (DAMASKSpectral.exe DAMASK_CONSTITUTIVE)

View File

@ -1,15 +0,0 @@
# group sources
set (KINEMATICS "kinematics_cleavage_opening"
"kinematics_slipplane_opening"
"kinematics_thermal_expansion"
"kinematics_vacancy_strain"
"kinematics_hydrogen_strain"
)
# compile module and cumulatively link the
# compiled libraries
foreach (p ${KINEMATICS})
add_library (${p} "${p}.f90")
target_link_libraries(${p} DAMASK_DRIVERS)
add_library (DAMASK_DRIVERS ALIAS ${p})
endforeach (p)

View File

@ -1,303 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Luv Sharma, Max-Planck-Institut fŸr Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut fŸr Eisenforschung GmbH
!> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_cleavage_opening
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
kinematics_cleavage_opening_sizePostResults, & !< cumulative size of post results
kinematics_cleavage_opening_offset, & !< which kinematics is my current damage mechanism?
kinematics_cleavage_opening_instance !< instance of damage kinematics mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
kinematics_cleavage_opening_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
kinematics_cleavage_opening_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
kinematics_cleavage_opening_Noutput !< number of outputs per instance of this damage
integer(pInt), dimension(:), allocatable, private :: &
kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems
integer(pInt), dimension(:,:), allocatable, private :: &
kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family
real(pReal), dimension(:), allocatable, private :: &
kinematics_cleavage_opening_sdot_0, &
kinematics_cleavage_opening_N
real(pReal), dimension(:,:), allocatable, private :: &
kinematics_cleavage_opening_critDisp, &
kinematics_cleavage_opening_critLoad
public :: &
kinematics_cleavage_opening_init, &
kinematics_cleavage_opening_LiAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_kinematics, &
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_cleavage_opening_label, &
KINEMATICS_cleavage_opening_ID, &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
worldrank
use lattice, only: &
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,kinematics
integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_kinematics == KINEMATICS_cleavage_opening_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_cleavage_opening_offset(material_Nphase), source=0_pInt)
allocate(kinematics_cleavage_opening_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
kinematics_cleavage_opening_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_cleavage_opening_ID)
do kinematics = 1, phase_Nkinematics(phase)
if (phase_kinematics(kinematics,phase) == kinematics_cleavage_opening_ID) &
kinematics_cleavage_opening_offset(phase) = kinematics
enddo
enddo
allocate(kinematics_cleavage_opening_sizePostResults(maxNinstance), source=0_pInt)
allocate(kinematics_cleavage_opening_sizePostResult(maxval(phase_Noutput),maxNinstance), source=0_pInt)
allocate(kinematics_cleavage_opening_output(maxval(phase_Noutput),maxNinstance))
kinematics_cleavage_opening_output = ''
allocate(kinematics_cleavage_opening_Noutput(maxNinstance), source=0_pInt)
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_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt)
allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0_pInt)
allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_cleavage_opening_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = kinematics_cleavage_opening_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('anisobrittle_sdot0')
kinematics_cleavage_opening_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('anisobrittle_ratesensitivity')
kinematics_cleavage_opening_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('ncleavage') !
Nchunks_CleavageFamilies = chunkPos(1) - 1_pInt
do j = 1_pInt, Nchunks_CleavageFamilies
kinematics_cleavage_opening_Ncleavage(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
enddo
case ('anisobrittle_criticaldisplacement')
do j = 1_pInt, Nchunks_CleavageFamilies
kinematics_cleavage_opening_critDisp(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
case ('anisobrittle_criticalload')
do j = 1_pInt, Nchunks_CleavageFamilies
kinematics_cleavage_opening_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
end select
endif; endif
enddo parsingFile
!--------------------------------------------------------------------------------------------------
! sanity checks
sanityChecks: do phase = 1_pInt, material_Nphase
myPhase: if (any(phase_kinematics(:,phase) == KINEMATICS_cleavage_opening_ID)) then
instance = kinematics_cleavage_opening_instance(phase)
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,phase),& ! limit active cleavage systems per family to min of available and requested
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance))
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) &
call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critDisp(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critLoad(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
endif myPhase
enddo sanityChecks
end subroutine kinematics_cleavage_opening_init
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el)
use prec, only: &
tol_math_check
use material, only: &
phaseAt, phasememberAt, &
material_homog, &
damage, &
damageMapping
use lattice, only: &
lattice_Scleavage, &
lattice_Scleavage_v, &
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar3333 !< derivative of Ld with respect to Tstar (4th-order tensor)
integer(pInt) :: &
phase, &
constituent, &
instance, &
homog, damageOffset, &
f, i, index_myFamily, k, l, m, n
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = kinematics_cleavage_opening_instance(phase)
homog = material_homog(ip,el)
damageOffset = damageMapping(homog)%p(ip,el)
Ld = 0.0_pReal
dLd_dTstar3333 = 0.0_pReal
do f = 1_pInt,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* &
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
udotd = &
sign(1.0_pReal,traction_d)* &
kinematics_cleavage_opening_sdot_0(instance)* &
(max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
if (abs(udotd) > tol_math_check) then
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)/ &
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) &
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,1,index_myFamily+i,phase)
endif
udott = &
sign(1.0_pReal,traction_t)* &
kinematics_cleavage_opening_sdot_0(instance)* &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
if (abs(udott) > tol_math_check) then
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)/ &
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) &
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,2,index_myFamily+i,phase)
endif
udotn = &
sign(1.0_pReal,traction_n)* &
kinematics_cleavage_opening_sdot_0(instance)* &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
if (abs(udotn) > tol_math_check) then
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)/ &
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) &
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,3,index_myFamily+i,phase)
endif
enddo
enddo
end subroutine kinematics_cleavage_opening_LiAndItsTangent
end module kinematics_cleavage_opening

View File

@ -1,264 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incorporating kinematics resulting from interstitial hydrogen
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_hydrogen_strain
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
kinematics_hydrogen_strain_sizePostResults, & !< cumulative size of post results
kinematics_hydrogen_strain_offset, & !< which kinematics is my current damage mechanism?
kinematics_hydrogen_strain_instance !< instance of damage kinematics mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
kinematics_hydrogen_strain_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
kinematics_hydrogen_strain_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
kinematics_hydrogen_strain_Noutput !< number of outputs per instance of this damage
real(pReal), dimension(:), allocatable, private :: &
kinematics_hydrogen_strain_coeff
public :: &
kinematics_hydrogen_strain_init, &
kinematics_hydrogen_strain_initialStrain, &
kinematics_hydrogen_strain_LiAndItsTangent, &
kinematics_hydrogen_strain_ChemPotAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_hydrogen_strain_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_kinematics, &
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_hydrogen_strain_label, &
KINEMATICS_hydrogen_strain_ID, &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
worldrank
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,kinematics
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_hydrogen_strain_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_kinematics == KINEMATICS_hydrogen_strain_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_hydrogen_strain_offset(material_Nphase), source=0_pInt)
allocate(kinematics_hydrogen_strain_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
kinematics_hydrogen_strain_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_hydrogen_strain_ID)
do kinematics = 1, phase_Nkinematics(phase)
if (phase_kinematics(kinematics,phase) == kinematics_hydrogen_strain_ID) &
kinematics_hydrogen_strain_offset(phase) = kinematics
enddo
enddo
allocate(kinematics_hydrogen_strain_sizePostResults(maxNinstance), source=0_pInt)
allocate(kinematics_hydrogen_strain_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(kinematics_hydrogen_strain_output(maxval(phase_Noutput),maxNinstance))
kinematics_hydrogen_strain_output = ''
allocate(kinematics_hydrogen_strain_Noutput(maxNinstance), source=0_pInt)
allocate(kinematics_hydrogen_strain_coeff(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_hydrogen_strain_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = kinematics_hydrogen_strain_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('hydrogen_strain_coeff')
kinematics_hydrogen_strain_coeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif; endif
enddo parsingFile
end subroutine kinematics_hydrogen_strain_init
!--------------------------------------------------------------------------------------------------
!> @brief report initial hydrogen strain based on current hydrogen conc deviation from
!> equillibrium (0)
!--------------------------------------------------------------------------------------------------
pure function kinematics_hydrogen_strain_initialStrain(ipc, ip, el)
use math, only: &
math_I3
use material, only: &
material_phase, &
material_homog, &
hydrogenConc, &
hydrogenfluxMapping
use lattice, only: &
lattice_equilibriumHydrogenConcentration
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
kinematics_hydrogen_strain_initialStrain !< initial thermal strain (should be small strain, though)
integer(pInt) :: &
phase, &
homog, offset, instance
phase = material_phase(ipc,ip,el)
instance = kinematics_hydrogen_strain_instance(phase)
homog = material_homog(ip,el)
offset = hydrogenfluxMapping(homog)%p(ip,el)
kinematics_hydrogen_strain_initialStrain = &
(hydrogenConc(homog)%p(offset) - lattice_equilibriumHydrogenConcentration(phase)) * &
kinematics_hydrogen_strain_coeff(instance)* math_I3
end function kinematics_hydrogen_strain_initialStrain
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_hydrogen_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el)
use material, only: &
material_phase, &
material_homog, &
hydrogenConc, &
hydrogenConcRate, &
hydrogenfluxMapping
use math, only: &
math_I3
use lattice, only: &
lattice_equilibriumHydrogenConcentration
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar3333 !< derivative of Li with respect to Tstar (4th-order tensor)
integer(pInt) :: &
phase, &
instance, &
homog, offset
real(pReal) :: &
Ch, ChEq, ChDot
phase = material_phase(ipc,ip,el)
instance = kinematics_hydrogen_strain_instance(phase)
homog = material_homog(ip,el)
offset = hydrogenfluxMapping(homog)%p(ip,el)
Ch = hydrogenConc(homog)%p(offset)
ChDot = hydrogenConcRate(homog)%p(offset)
ChEq = lattice_equilibriumHydrogenConcentration(phase)
Li = ChDot*math_I3* &
kinematics_hydrogen_strain_coeff(instance)/ &
(1.0_pReal + kinematics_hydrogen_strain_coeff(instance)*(Ch - ChEq))
dLi_dTstar3333 = 0.0_pReal
end subroutine kinematics_hydrogen_strain_LiAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief contains the kinematic contribution to hydrogen chemical potential
!--------------------------------------------------------------------------------------------------
subroutine kinematics_hydrogen_strain_ChemPotAndItsTangent(ChemPot, dChemPot_dCh, Tstar_v, Fi0, Fi, ipc, ip, el)
use material, only: &
material_phase
use math, only: &
math_inv33, &
math_mul33x33, &
math_Mandel6to33, &
math_transpose33
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(6) :: &
Tstar_v
real(pReal), intent(in), dimension(3,3) :: &
Fi0, Fi
real(pReal), intent(out) :: &
ChemPot, dChemPot_dCh
integer(pInt) :: &
phase, &
instance
phase = material_phase(ipc,ip,el)
instance = kinematics_hydrogen_strain_instance(phase)
ChemPot = -kinematics_hydrogen_strain_coeff(instance)* &
sum(math_mul33x33(Fi,math_Mandel6to33(Tstar_v))* &
math_mul33x33(math_mul33x33(Fi,math_inv33(Fi0)),Fi))
dChemPot_dCh = 0.0_pReal
end subroutine kinematics_hydrogen_strain_ChemPotAndItsTangent
end module kinematics_hydrogen_strain

View File

@ -1,323 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incorporating kinematics resulting from opening of slip planes
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_slipplane_opening
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
kinematics_slipplane_opening_sizePostResults, & !< cumulative size of post results
kinematics_slipplane_opening_offset, & !< which kinematics is my current damage mechanism?
kinematics_slipplane_opening_instance !< instance of damage kinematics mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
kinematics_slipplane_opening_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
kinematics_slipplane_opening_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
kinematics_slipplane_opening_Noutput !< number of outputs per instance of this damage
integer(pInt), dimension(:), allocatable, private :: &
kinematics_slipplane_opening_totalNslip !< total number of slip systems
integer(pInt), dimension(:,:), allocatable, private :: &
kinematics_slipplane_opening_Nslip !< number of slip systems per family
real(pReal), dimension(:), allocatable, private :: &
kinematics_slipplane_opening_sdot_0, &
kinematics_slipplane_opening_N
real(pReal), dimension(:,:), allocatable, private :: &
kinematics_slipplane_opening_critPlasticStrain, &
kinematics_slipplane_opening_critLoad
public :: &
kinematics_slipplane_opening_init, &
kinematics_slipplane_opening_LiAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_kinematics, &
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_slipplane_opening_label, &
KINEMATICS_slipplane_opening_ID, &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
worldrank
use lattice, only: &
lattice_maxNslipFamily, &
lattice_NslipSystem
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,kinematics
integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_kinematics == KINEMATICS_slipplane_opening_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_slipplane_opening_offset(material_Nphase), source=0_pInt)
allocate(kinematics_slipplane_opening_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
kinematics_slipplane_opening_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_slipplane_opening_ID)
do kinematics = 1, phase_Nkinematics(phase)
if (phase_kinematics(kinematics,phase) == kinematics_slipplane_opening_ID) &
kinematics_slipplane_opening_offset(phase) = kinematics
enddo
enddo
allocate(kinematics_slipplane_opening_sizePostResults(maxNinstance), source=0_pInt)
allocate(kinematics_slipplane_opening_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(kinematics_slipplane_opening_output(maxval(phase_Noutput),maxNinstance))
kinematics_slipplane_opening_output = ''
allocate(kinematics_slipplane_opening_Noutput(maxNinstance), source=0_pInt)
allocate(kinematics_slipplane_opening_critLoad(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal)
allocate(kinematics_slipplane_opening_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal)
allocate(kinematics_slipplane_opening_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt)
allocate(kinematics_slipplane_opening_totalNslip(maxNinstance), source=0_pInt)
allocate(kinematics_slipplane_opening_N(maxNinstance), source=0.0_pReal)
allocate(kinematics_slipplane_opening_sdot_0(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_slipplane_opening_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = kinematics_slipplane_opening_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('nslip') !
Nchunks_SlipFamilies = chunkPos(1) - 1_pInt
do j = 1_pInt, Nchunks_SlipFamilies
kinematics_slipplane_opening_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
enddo
case ('anisoductile_sdot0')
kinematics_slipplane_opening_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('anisoductile_criticalplasticstrain')
do j = 1_pInt, Nchunks_SlipFamilies
kinematics_slipplane_opening_critPlasticStrain(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
case ('anisoductile_ratesensitivity')
kinematics_slipplane_opening_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('anisoductile_criticalload')
do j = 1_pInt, Nchunks_SlipFamilies
kinematics_slipplane_opening_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
end select
endif; endif
enddo parsingFile
!--------------------------------------------------------------------------------------------------
! sanity checks
sanityChecks: do phase = 1_pInt, material_Nphase
myPhase: if (any(phase_kinematics(:,phase) == KINEMATICS_slipplane_opening_ID)) then
instance = kinematics_slipplane_opening_instance(phase)
kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance) = &
min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active cleavage systems per family to min of available and requested
kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance))
kinematics_slipplane_opening_totalNslip(instance) = sum(kinematics_slipplane_opening_Nslip(:,instance))
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//')')
if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')')
if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')')
endif myPhase
enddo sanityChecks
end subroutine kinematics_slipplane_opening_init
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el)
use prec, only: &
tol_math_check
use lattice, only: &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_sd, &
lattice_st, &
lattice_sn
use material, only: &
phaseAt, phasememberAt, &
material_homog, &
damage, &
damageMapping
use math, only: &
math_Plain3333to99, &
math_I3, &
math_identity4th, &
math_symmetric33, &
math_Mandel33to6, &
math_tensorproduct33, &
math_det33, &
math_mul33x33
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar3333 !< derivative of Ld with respect to Tstar (4th-order tensor)
real(pReal), dimension(3,3) :: &
projection_d, projection_t, projection_n !< projection modes 3x3 tensor
real(pReal), dimension(6) :: &
projection_d_v, projection_t_v, projection_n_v !< projection modes 3x3 vector
integer(pInt) :: &
phase, &
constituent, &
instance, &
homog, damageOffset, &
f, i, index_myFamily, k, l, m, n
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = kinematics_slipplane_opening_instance(phase)
homog = material_homog(ip,el)
damageOffset = damageMapping(homog)%p(ip,el)
Ld = 0.0_pReal
dLd_dTstar3333 = 0.0_pReal
do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,kinematics_slipplane_opening_Nslip(f,instance) ! process each (active) slip system in family
projection_d = math_tensorproduct33(lattice_sd(1:3,index_myFamily+i,phase),&
lattice_sn(1:3,index_myFamily+i,phase))
projection_t = math_tensorproduct33(lattice_st(1:3,index_myFamily+i,phase),&
lattice_sn(1:3,index_myFamily+i,phase))
projection_n = math_tensorproduct33(lattice_sn(1:3,index_myFamily+i,phase),&
lattice_sn(1:3,index_myFamily+i,phase))
projection_d_v(1:6) = math_Mandel33to6(math_symmetric33(projection_d(1:3,1:3)))
projection_t_v(1:6) = math_Mandel33to6(math_symmetric33(projection_t(1:3,1:3)))
projection_n_v(1:6) = math_Mandel33to6(math_symmetric33(projection_n(1:3,1:3)))
traction_d = dot_product(Tstar_v,projection_d_v(1:6))
traction_t = dot_product(Tstar_v,projection_t_v(1:6))
traction_n = dot_product(Tstar_v,projection_n_v(1:6))
traction_crit = kinematics_slipplane_opening_critLoad(f,instance)* &
damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage
udotd = &
sign(1.0_pReal,traction_d)* &
kinematics_slipplane_opening_sdot_0(instance)* &
(abs(traction_d)/traction_crit - &
abs(traction_d)/kinematics_slipplane_opening_critLoad(f,instance))**kinematics_slipplane_opening_N(instance)
if (abs(udotd) > tol_math_check) then
Ld = Ld + udotd*projection_d
dudotd_dt = udotd*kinematics_slipplane_opening_N(instance)/traction_d
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudotd_dt*projection_d(k,l)*projection_d(m,n)
endif
udott = &
sign(1.0_pReal,traction_t)* &
kinematics_slipplane_opening_sdot_0(instance)* &
(abs(traction_t)/traction_crit - &
abs(traction_t)/kinematics_slipplane_opening_critLoad(f,instance))**kinematics_slipplane_opening_N(instance)
if (abs(udott) > tol_math_check) then
Ld = Ld + udott*projection_t
dudott_dt = udott*kinematics_slipplane_opening_N(instance)/traction_t
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudott_dt*projection_t(k,l)*projection_t(m,n)
endif
udotn = &
kinematics_slipplane_opening_sdot_0(instance)* &
(max(0.0_pReal,traction_n)/traction_crit - &
max(0.0_pReal,traction_n)/kinematics_slipplane_opening_critLoad(f,instance))**kinematics_slipplane_opening_N(instance)
if (abs(udotn) > tol_math_check) then
Ld = Ld + udotn*projection_n
dudotn_dt = udotn*kinematics_slipplane_opening_N(instance)/traction_n
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
dudotn_dt*projection_n(k,l)*projection_n(m,n)
endif
enddo
enddo
end subroutine kinematics_slipplane_opening_LiAndItsTangent
end module kinematics_slipplane_opening

View File

@ -1,228 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incorporating kinematics resulting from thermal expansion
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_thermal_expansion
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
kinematics_thermal_expansion_sizePostResults, & !< cumulative size of post results
kinematics_thermal_expansion_offset, & !< which kinematics is my current damage mechanism?
kinematics_thermal_expansion_instance !< instance of damage kinematics mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
kinematics_thermal_expansion_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
kinematics_thermal_expansion_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
kinematics_thermal_expansion_Noutput !< number of outputs per instance of this damage
! enum, bind(c) ! ToDo kinematics need state machinery to deal with sizePostResult
! enumerator :: undefined_ID, & ! possible remedy is to decouple having state vars from having output
! thermalexpansionrate_ID ! which means to separate user-defined types tState + tOutput...
! end enum
public :: &
kinematics_thermal_expansion_init, &
kinematics_thermal_expansion_initialStrain, &
kinematics_thermal_expansion_LiAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_kinematics, &
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_thermal_expansion_label, &
KINEMATICS_thermal_expansion_ID, &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
worldrank
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,kinematics
character(len=65536) :: &
tag = '', &
output = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_thermal_expansion_offset(material_Nphase), source=0_pInt)
allocate(kinematics_thermal_expansion_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
kinematics_thermal_expansion_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_thermal_expansion_ID)
do kinematics = 1, phase_Nkinematics(phase)
if (phase_kinematics(kinematics,phase) == kinematics_thermal_expansion_ID) &
kinematics_thermal_expansion_offset(phase) = kinematics
enddo
enddo
allocate(kinematics_thermal_expansion_sizePostResults(maxNinstance), source=0_pInt)
allocate(kinematics_thermal_expansion_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(kinematics_thermal_expansion_output(maxval(phase_Noutput),maxNinstance))
kinematics_thermal_expansion_output = ''
allocate(kinematics_thermal_expansion_Noutput(maxNinstance), source=0_pInt)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_thermal_expansion_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = kinematics_thermal_expansion_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key...
select case(tag)
! case ('(output)')
! output = IO_lc(IO_stringValue(line,chunkPos,2_pInt)) ! ...and corresponding output
! select case(output)
! case ('thermalexpansionrate')
! kinematics_thermal_expansion_Noutput(instance) = kinematics_thermal_expansion_Noutput(instance) + 1_pInt
! kinematics_thermal_expansion_outputID(kinematics_thermal_expansion_Noutput(instance),instance) = &
! thermalexpansionrate_ID
! kinematics_thermal_expansion_output(kinematics_thermal_expansion_Noutput(instance),instance) = output
! ToDo add sizePostResult loop afterwards...
end select
endif; endif
enddo parsingFile
end subroutine kinematics_thermal_expansion_init
!--------------------------------------------------------------------------------------------------
!> @brief report initial thermal strain based on current temperature deviation from reference
!--------------------------------------------------------------------------------------------------
pure function kinematics_thermal_expansion_initialStrain(ipc, ip, el)
use material, only: &
material_phase, &
material_homog, &
temperature, &
thermalMapping
use lattice, only: &
lattice_thermalExpansion33, &
lattice_referenceTemperature
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though)
integer(pInt) :: &
phase, &
homog, offset
phase = material_phase(ipc,ip,el)
homog = material_homog(ip,el)
offset = thermalMapping(homog)%p(ip,el)
kinematics_thermal_expansion_initialStrain = &
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase)) * &
lattice_thermalExpansion33(1:3,1:3,phase)
end function kinematics_thermal_expansion_initialStrain
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el)
use material, only: &
material_phase, &
material_homog, &
temperature, &
temperatureRate, &
thermalMapping
use lattice, only: &
lattice_thermalExpansion33, &
lattice_referenceTemperature
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar3333 !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
integer(pInt) :: &
phase, &
homog, offset
real(pReal) :: &
T, TRef, TDot
phase = material_phase(ipc,ip,el)
homog = material_homog(ip,el)
offset = thermalMapping(homog)%p(ip,el)
T = temperature(homog)%p(offset)
TDot = temperatureRate(homog)%p(offset)
TRef = lattice_referenceTemperature(phase)
Li = TDot* &
lattice_thermalExpansion33(1:3,1:3,phase)/ &
(1.0_pReal + lattice_thermalExpansion33(1:3,1:3,phase)*(T - TRef))
dLi_dTstar3333 = 0.0_pReal
end subroutine kinematics_thermal_expansion_LiAndItsTangent
end module kinematics_thermal_expansion

View File

@ -1,265 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incorporating kinematics resulting from vacancy point defects
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_vacancy_strain
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
kinematics_vacancy_strain_sizePostResults, & !< cumulative size of post results
kinematics_vacancy_strain_offset, & !< which kinematics is my current damage mechanism?
kinematics_vacancy_strain_instance !< instance of damage kinematics mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
kinematics_vacancy_strain_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
kinematics_vacancy_strain_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
kinematics_vacancy_strain_Noutput !< number of outputs per instance of this damage
real(pReal), dimension(:), allocatable, private :: &
kinematics_vacancy_strain_coeff
public :: &
kinematics_vacancy_strain_init, &
kinematics_vacancy_strain_initialStrain, &
kinematics_vacancy_strain_LiAndItsTangent, &
kinematics_vacancy_strain_ChemPotAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_vacancy_strain_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_kinematics, &
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_vacancy_strain_label, &
KINEMATICS_vacancy_strain_ID, &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
worldrank
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,kinematics
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_vacancy_strain_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_kinematics == KINEMATICS_vacancy_strain_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_vacancy_strain_offset(material_Nphase), source=0_pInt)
allocate(kinematics_vacancy_strain_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
kinematics_vacancy_strain_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_vacancy_strain_ID)
do kinematics = 1, phase_Nkinematics(phase)
if (phase_kinematics(kinematics,phase) == kinematics_vacancy_strain_ID) &
kinematics_vacancy_strain_offset(phase) = kinematics
enddo
enddo
allocate(kinematics_vacancy_strain_sizePostResults(maxNinstance), source=0_pInt)
allocate(kinematics_vacancy_strain_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(kinematics_vacancy_strain_output(maxval(phase_Noutput),maxNinstance))
kinematics_vacancy_strain_output = ''
allocate(kinematics_vacancy_strain_Noutput(maxNinstance), source=0_pInt)
allocate(kinematics_vacancy_strain_coeff(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_vacancy_strain_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = kinematics_vacancy_strain_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('vacancy_strain_coeff')
kinematics_vacancy_strain_coeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif; endif
enddo parsingFile
end subroutine kinematics_vacancy_strain_init
!--------------------------------------------------------------------------------------------------
!> @brief report initial vacancy strain based on current vacancy conc deviation from equillibrium
!--------------------------------------------------------------------------------------------------
pure function kinematics_vacancy_strain_initialStrain(ipc, ip, el)
use math, only: &
math_I3
use material, only: &
material_phase, &
material_homog, &
vacancyConc, &
vacancyfluxMapping
use lattice, only: &
lattice_equilibriumVacancyConcentration
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
kinematics_vacancy_strain_initialStrain !< initial thermal strain (should be small strain, though)
integer(pInt) :: &
phase, &
homog, offset, instance
phase = material_phase(ipc,ip,el)
instance = kinematics_vacancy_strain_instance(phase)
homog = material_homog(ip,el)
offset = vacancyfluxMapping(homog)%p(ip,el)
kinematics_vacancy_strain_initialStrain = &
(vacancyConc(homog)%p(offset) - lattice_equilibriumVacancyConcentration(phase)) * &
kinematics_vacancy_strain_coeff(instance)* math_I3
end function kinematics_vacancy_strain_initialStrain
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_vacancy_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el)
use material, only: &
material_phase, &
material_homog, &
vacancyConc, &
vacancyConcRate, &
vacancyfluxMapping
use math, only: &
math_I3
use lattice, only: &
lattice_equilibriumVacancyConcentration
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar3333 !< derivative of Li with respect to Tstar (4th-order tensor)
integer(pInt) :: &
phase, &
instance, &
homog, offset
real(pReal) :: &
Cv, CvEq, CvDot
phase = material_phase(ipc,ip,el)
instance = kinematics_vacancy_strain_instance(phase)
homog = material_homog(ip,el)
offset = vacancyfluxMapping(homog)%p(ip,el)
Cv = vacancyConc(homog)%p(offset)
CvDot = vacancyConcRate(homog)%p(offset)
CvEq = lattice_equilibriumvacancyConcentration(phase)
Li = CvDot*math_I3* &
kinematics_vacancy_strain_coeff(instance)/ &
(1.0_pReal + kinematics_vacancy_strain_coeff(instance)*(Cv - CvEq))
dLi_dTstar3333 = 0.0_pReal
end subroutine kinematics_vacancy_strain_LiAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief contains the kinematic contribution to vacancy chemical potential
!--------------------------------------------------------------------------------------------------
subroutine kinematics_vacancy_strain_ChemPotAndItsTangent(ChemPot, dChemPot_dCv, Tstar_v, Fi0, Fi, ipc, ip, el)
use material, only: &
material_phase
use math, only: &
math_inv33, &
math_mul33x33, &
math_Mandel6to33, &
math_transpose33
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(6) :: &
Tstar_v
real(pReal), intent(in), dimension(3,3) :: &
Fi0, Fi
real(pReal), intent(out) :: &
ChemPot, dChemPot_dCv
integer(pInt) :: &
phase, &
instance
phase = material_phase(ipc,ip,el)
instance = kinematics_vacancy_strain_instance(phase)
ChemPot = -kinematics_vacancy_strain_coeff(instance)* &
sum(math_mul33x33(Fi,math_Mandel6to33(Tstar_v))* &
math_mul33x33(math_mul33x33(Fi,math_inv33(Fi0)),Fi))
dChemPot_dCv = 0.0_pReal
end subroutine kinematics_vacancy_strain_ChemPotAndItsTangent
end module kinematics_vacancy_strain

View File

@ -1,29 +0,0 @@
# group sources
set (PLASTIC "plastic_dislotwin"
"plastic_disloUCLA"
"plastic_isotropic"
"plastic_j2"
"plastic_phenopowerlaw"
"plastic_titanmod"
"plastic_nonlocal"
"plastic_none"
"plastic_phenoplus"
)
# compile module and cumulatively link the
# compiled libraries
add_library (DAMASK_PLASTIC "plastic_dislotwin.f90"
"plastic_disloUCLA.f90"
"plastic_isotropic.f90"
"plastic_j2.f90"
"plastic_phenopowerlaw.f90"
"plastic_titanmod.f90"
"plastic_nonlocal.f90"
"plastic_none.f90"
"plastic_phenoplus.f90")
target_link_libraries(DAMASK_PLASTIC DAMASK_DRIVERS)
# foreach (p ${PLASTIC})
# add_library (${p} "${p}.f90")
# target_link_libraries(${p} DAMASK_DRIVERS)
# add_library (DAMASK_DRIVERS ALIAS ${p})
# endforeach (p)

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -1,678 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for isotropic (ISOTROPIC) plasticity
!> @details Isotropic (ISOTROPIC) Plasticity which resembles the phenopowerlaw plasticity without
!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
module plastic_isotropic
#ifdef HDF
use hdf5, only: &
HID_T
#endif
use prec, only: &
pReal,&
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_isotropic_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: &
plastic_isotropic_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_isotropic_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
plastic_isotropic_Noutput !< number of outputs per instance
enum, bind(c)
enumerator :: undefined_ID, &
flowstress_ID, &
strainrate_ID
end enum
type, private :: tParameters !< container type for internal constitutive parameters
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
real(pReal) :: &
fTaylor, &
tau0, &
gdot0, &
n, &
h0, &
h0_slopeLnRate, &
tausat, &
a, &
aTolFlowstress, &
aTolShear , &
tausat_SinhFitA, &
tausat_SinhFitB, &
tausat_SinhFitC, &
tausat_SinhFitD
logical :: &
dilatation
end type
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type, private :: tIsotropicState !< internal state aliases
real(pReal), pointer, dimension(:) :: & ! scalars along NipcMyInstance
flowstress, &
accumulatedShear
end type
type, private :: tIsotropicAbsTol !< internal alias for abs tolerance in state
real(pReal), pointer :: & ! scalars along NipcMyInstance
flowstress, &
accumulatedShear
end type
type(tIsotropicState), allocatable, dimension(:), private :: & !< state aliases per instance
state, &
state0, &
dotState
type(tIsotropicAbsTol), allocatable, dimension(:), private :: & !< state aliases per instance
stateAbsTol
public :: &
plastic_isotropic_init, &
plastic_isotropic_LpAndItsTangent, &
plastic_isotropic_LiAndItsTangent, &
plastic_isotropic_dotState, &
plastic_isotropic_postResults
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelBasic
use numerics, only: &
analyticJaco, &
worldrank, &
numerics_integrator
use math, only: &
math_Mandel3333to66, &
math_Voigt66to3333
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_plasticity, &
phase_plasticityInstance, &
phase_Noutput, &
PLASTICITY_ISOTROPIC_label, &
PLASTICITY_ISOTROPIC_ID, &
material_phase, &
plasticState, &
MATERIAL_partPhase
use lattice
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: &
o, &
phase, &
instance, &
maxNinstance, &
mySize, &
sizeDotState, &
sizeState, &
sizeDeltaState
character(len=65536) :: &
tag = '', &
outputtag = '', &
line = '', &
extmsg = ''
integer(pInt) :: NipcMyPhase
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_ISOTROPIC_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(plastic_isotropic_sizePostResults(maxNinstance), source=0_pInt)
allocate(plastic_isotropic_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
allocate(plastic_isotropic_output(maxval(phase_Noutput), maxNinstance))
plastic_isotropic_output = ''
allocate(plastic_isotropic_Noutput(maxNinstance), source=0_pInt)
allocate(param(maxNinstance)) ! one container of parameters per instance
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
phase = phase + 1_pInt ! advance section counter
if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then
instance = phase_plasticityInstance(phase)
endif
cycle ! skip to next line
endif
if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
allocate(param(instance)%outputID(phase_Noutput(phase))) ! allocate space for IDs of every requested output
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
extmsg = trim(tag)//' ('//PLASTICITY_ISOTROPIC_label//')' ! prepare error message identifier
select case(tag)
case ('(output)')
outputtag = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
select case(outputtag)
case ('flowstress')
plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
param(instance)%outputID (plastic_isotropic_Noutput(instance)) = flowstress_ID
plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag
case ('strainrate')
plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
param(instance)%outputID (plastic_isotropic_Noutput(instance)) = strainrate_ID
plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag
end select
case ('/dilatation/')
param(instance)%dilatation = .true.
case ('tau0')
param(instance)%tau0 = IO_floatValue(line,chunkPos,2_pInt)
if (param(instance)%tau0 < 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg)
case ('gdot0')
param(instance)%gdot0 = IO_floatValue(line,chunkPos,2_pInt)
if (param(instance)%gdot0 <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg)
case ('n')
param(instance)%n = IO_floatValue(line,chunkPos,2_pInt)
if (param(instance)%n <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg)
case ('h0')
param(instance)%h0 = IO_floatValue(line,chunkPos,2_pInt)
case ('h0_slope','slopelnrate')
param(instance)%h0_slopeLnRate = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat')
param(instance)%tausat = IO_floatValue(line,chunkPos,2_pInt)
if (param(instance)%tausat <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg)
case ('tausat_sinhfita')
param(instance)%tausat_SinhFitA = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitb')
param(instance)%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitc')
param(instance)%tausat_SinhFitC = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitd')
param(instance)%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt)
case ('a', 'w0')
param(instance)%a = IO_floatValue(line,chunkPos,2_pInt)
if (param(instance)%a <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg)
case ('taylorfactor')
param(instance)%fTaylor = IO_floatValue(line,chunkPos,2_pInt)
if (param(instance)%fTaylor <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg)
case ('atol_flowstress')
param(instance)%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt)
if (param(instance)%aTolFlowstress <= 0.0_pReal) call IO_error(211_pInt,ext_msg=extmsg)
case ('atol_shear')
param(instance)%aTolShear = IO_floatValue(line,chunkPos,2_pInt)
case default
end select
endif; endif
enddo parsingFile
allocate(state(maxNinstance)) ! internal state aliases
allocate(state0(maxNinstance))
allocate(dotState(maxNinstance))
allocate(stateAbsTol(maxNinstance))
initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop over every plasticity
myPhase: if (phase_plasticity(phase) == PLASTICITY_isotropic_ID) then ! isolate instances of own constitutive description
NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc)
instance = phase_plasticityInstance(phase)
!--------------------------------------------------------------------------------------------------
! sanity checks
if (param(instance)%aTolShear <= 0.0_pReal) &
param(instance)%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
select case(param(instance)%outputID(o))
case(flowstress_ID,strainrate_ID)
mySize = 1_pInt
case default
end select
outputFound: if (mySize > 0_pInt) then
plastic_isotropic_sizePostResult(o,instance) = mySize
plastic_isotropic_sizePostResults(instance) = &
plastic_isotropic_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
!--------------------------------------------------------------------------------------------------
! allocate state arrays
sizeState = 2_pInt ! flowstress, accumulated_shear
sizeDotState = sizeState ! both evolve
sizeDeltaState = 0_pInt ! no sudden jumps in state
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%sizePostResults = plastic_isotropic_sizePostResults(instance)
plasticState(phase)%nSlip = 1
plasticState(phase)%nTwin = 0
plasticState(phase)%nTrans= 0
allocate(plasticState(phase)%aTolState ( sizeState))
allocate(plasticState(phase)%state0 ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%partionedState0 ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%subState0 ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%state ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup ( sizeState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NipcMyPhase),source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal)
!--------------------------------------------------------------------------------------------------
! globally required state aliases
plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase)
plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases
state(instance)%flowstress => plasticState(phase)%state (1,1:NipcMyPhase)
state0(instance)%flowstress => plasticState(phase)%state0 (1,1:NipcMyPhase)
dotState(instance)%flowstress => plasticState(phase)%dotState (1,1:NipcMyPhase)
stateAbsTol(instance)%flowstress => plasticState(phase)%aTolState(1)
state(instance)%accumulatedShear => plasticState(phase)%state (2,1:NipcMyPhase)
state0(instance)%accumulatedShear => plasticState(phase)%state0 (2,1:NipcMyPhase)
dotState(instance)%accumulatedShear => plasticState(phase)%dotState (2,1:NipcMyPhase)
stateAbsTol(instance)%accumulatedShear => plasticState(phase)%aTolState(2)
!--------------------------------------------------------------------------------------------------
! init state
state0(instance)%flowstress = param(instance)%tau0
state0(instance)%accumulatedShear = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! init absolute state tolerances
stateAbsTol(instance)%flowstress = param(instance)%aTolFlowstress
stateAbsTol(instance)%accumulatedShear = param(instance)%aTolShear
endif myPhase
enddo initializeInstances
end subroutine plastic_isotropic_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective, &
debug_e, &
debug_i, &
debug_g
use math, only: &
math_mul6x6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_deviatoric33, &
math_mul33xx33, &
math_transpose33
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar_3333 !< derivative of Lp with respect to Tstar as 4th order tensor
real(pReal) :: &
gamma_dot, & !< strainrate
norm_Tstar_dev, & !< euclidean norm of Tstar_dev
squarenorm_Tstar_dev !< square of the euclidean norm of Tstar_dev
integer(pInt) :: &
instance, of, &
k, l, m, n
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !!
Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
norm_Tstar_dev = sqrt(squarenorm_Tstar_dev)
if (norm_Tstar_dev <= 0.0_pReal) then ! Tstar == 0 --> both Lp and dLp_dTstar are zero
Lp = 0.0_pReal
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = param(instance)%gdot0 &
* ( sqrt(1.5_pReal) * norm_Tstar_dev / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/param(instance)%fTaylor
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CONST isotropic >> at el ip g ',el,ip,ipc
write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CONST isotropic >> Tstar (dev) / MPa', &
math_transpose33(Tstar_dev_33(1:3,1:3))*1.0e-6_pReal
write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> norm Tstar / MPa', norm_Tstar_dev*1.0e-6_pReal
write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', gamma_dot
end if
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * &
Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
dLp_dTstar99 = math_Plain3333to99(gamma_dot / param(instance)%fTaylor * &
dLp_dTstar_3333 / norm_Tstar_dev)
end if
end subroutine plastic_isotropic_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_spherical33, &
math_mul33xx33
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
Li !< plastic velocity gradient
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(3,3) :: &
Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dTstar_3333 !< derivative of Li with respect to Tstar as 4th order tensor
real(pReal) :: &
gamma_dot, & !< strainrate
norm_Tstar_sph, & !< euclidean norm of Tstar_sph
squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph
integer(pInt) :: &
instance, of, &
k, l, m, n
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !!
Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33)
norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
if (param(instance)%dilatation) then
if (norm_Tstar_sph <= 0.0_pReal) then ! Tstar == 0 --> both Li and dLi_dTstar are zero
Li = 0.0_pReal
dLi_dTstar_3333 = 0.0_pReal
else
gamma_dot = param(instance)%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n
Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Li
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * &
Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * &
dLi_dTstar_3333 / norm_Tstar_sph
endif
endif
end subroutine plastic_isotropic_LiAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
gamma_dot, & !< strainrate
hardening, & !< hardening coefficient
saturation, & !< saturation flowstress
norm_Tstar_v !< euclidean norm of Tstar_dev
integer(pInt) :: &
instance, & !< instance of my instance (unique number of my constitutive model)
of !< shortcut notation for offset position in state array
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !!
!--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress
if (param(instance)%dilatation) then
norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
else
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
Tstar_dev_v(4:6) = Tstar_v(4:6)
norm_Tstar_v = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
end if
!--------------------------------------------------------------------------------------------------
! strain rate
gamma_dot = param(instance)%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
/ &!-----------------------------------------------------------------------------------
(param(instance)%fTaylor*state(instance)%flowstress(of) ))**param(instance)%n
!--------------------------------------------------------------------------------------------------
! hardening coefficient
if (abs(gamma_dot) > 1e-12_pReal) then
if (abs(param(instance)%tausat_SinhFitA) <= tiny(0.0_pReal)) then
saturation = param(instance)%tausat
else
saturation = ( param(instance)%tausat &
+ ( log( ( gamma_dot / param(instance)%tausat_SinhFitA&
)**(1.0_pReal / param(instance)%tausat_SinhFitD)&
+ sqrt( ( gamma_dot / param(instance)%tausat_SinhFitA &
)**(2.0_pReal / param(instance)%tausat_SinhFitD) &
+ 1.0_pReal ) &
) & ! asinh(K) = ln(K + sqrt(K^2 +1))
)**(1.0_pReal / param(instance)%tausat_SinhFitC) &
/ ( param(instance)%tausat_SinhFitB &
* (gamma_dot / param(instance)%gdot0)**(1.0_pReal / param(instance)%n) &
) &
)
endif
hardening = ( param(instance)%h0 + param(instance)%h0_slopeLnRate * log(gamma_dot) ) &
* abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**param(instance)%a &
* sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation)
else
hardening = 0.0_pReal
endif
dotState(instance)%flowstress (of) = hardening * gamma_dot
dotState(instance)%accumulatedShear(of) = gamma_dot
end subroutine plastic_isotropic_dotState
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use material, only: &
material_phase, &
plasticState, &
phaseAt, phasememberAt, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
plastic_isotropic_postResults
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
norm_Tstar_v ! euclidean norm of Tstar_dev
integer(pInt) :: &
instance, & !< instance of my instance (unique number of my constitutive model)
of, & !< shortcut notation for offset position in state array
c, &
o
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !!
!--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress
if (param(instance)%dilatation) then
norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
else
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
Tstar_dev_v(4:6) = Tstar_v(4:6)
norm_Tstar_v = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
end if
c = 0_pInt
plastic_isotropic_postResults = 0.0_pReal
outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
select case(param(instance)%outputID(o))
case (flowstress_ID)
plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of)
c = c + 1_pInt
case (strainrate_ID)
plastic_isotropic_postResults(c+1_pInt) = &
param(instance)%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
/ &!----------------------------------------------------------------------------------
(param(instance)%fTaylor * state(instance)%flowstress(of)) ) ** param(instance)%n
c = c + 1_pInt
end select
enddo outputsLoop
end function plastic_isotropic_postResults
end module plastic_isotropic

View File

@ -1,579 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for isotropic (J2) plasticity
!> @details Isotropic (J2) Plasticity which resembles the phenopowerlaw plasticity without
!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
module plastic_j2
#ifdef HDF
use hdf5, only: &
HID_T
#endif
use prec, only: &
pReal,&
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_j2_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: &
plastic_j2_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
plastic_j2_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
plastic_j2_Noutput !< number of outputs per instance
real(pReal), dimension(:), allocatable, private :: &
plastic_j2_fTaylor, & !< Taylor factor
plastic_j2_tau0, & !< initial plastic stress
plastic_j2_gdot0, & !< reference velocity
plastic_j2_n, & !< Visco-plastic parameter
!--------------------------------------------------------------------------------------------------
! h0 as function of h0 = A + B log (gammadot)
plastic_j2_h0, &
plastic_j2_h0_slopeLnRate, &
plastic_j2_tausat, & !< final plastic stress
plastic_j2_a, &
plastic_j2_aTolResistance, &
plastic_j2_aTolShear, &
!--------------------------------------------------------------------------------------------------
! tausat += (asinh((gammadot / SinhFitA)**(1 / SinhFitD)))**(1 / SinhFitC) / (SinhFitB * (gammadot / gammadot0)**(1/n))
plastic_j2_tausat_SinhFitA, & !< fitting parameter for normalized strain rate vs. stress function
plastic_j2_tausat_SinhFitB, & !< fitting parameter for normalized strain rate vs. stress function
plastic_j2_tausat_SinhFitC, & !< fitting parameter for normalized strain rate vs. stress function
plastic_j2_tausat_SinhFitD !< fitting parameter for normalized strain rate vs. stress function
enum, bind(c)
enumerator :: undefined_ID, &
flowstress_ID, &
strainrate_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
plastic_j2_outputID !< ID of each post result output
#ifdef HDF
type plastic_j2_tOutput
real(pReal), dimension(:), allocatable, private :: &
flowstress, &
strainrate
logical :: flowstressActive = .false., strainrateActive = .false. ! if we can write the output block wise, this is not needed anymore because we can do an if(allocated(xxx))
end type plastic_j2_tOutput
type(plastic_j2_tOutput), allocatable, dimension(:) :: plastic_j2_Output2
integer(HID_T), allocatable, dimension(:) :: outID
#endif
public :: &
plastic_j2_init, &
plastic_j2_LpAndItsTangent, &
plastic_j2_dotState, &
plastic_j2_postResults
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
#ifdef HDF
use hdf5
#endif
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelBasic
use numerics, only: &
analyticJaco, &
worldrank, &
numerics_integrator
use math, only: &
math_Mandel3333to66, &
math_Voigt66to3333
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_error, &
IO_timeStamp, &
#ifdef HDF
tempResults, &
HDF5_addGroup, &
HDF5_addScalarDataset,&
#endif
IO_EOF
use material, only: &
phase_plasticity, &
phase_plasticityInstance, &
phase_Noutput, &
PLASTICITY_J2_label, &
PLASTICITY_J2_ID, &
material_phase, &
plasticState, &
MATERIAL_partPhase
use lattice
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: &
o, &
phase, &
maxNinstance, &
instance, &
mySize, &
sizeDotState, &
sizeState, &
sizeDeltaState
character(len=65536) :: &
tag = '', &
line = ''
integer(pInt) :: NofMyPhase
#ifdef HDF
character(len=5) :: &
str1
integer(HID_T) :: ID,ID2,ID4
#endif
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_J2_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_J2_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#ifdef HDF
allocate(plastic_j2_Output2(maxNinstance))
allocate(outID(maxNinstance))
#endif
allocate(plastic_j2_sizePostResults(maxNinstance), source=0_pInt)
allocate(plastic_j2_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
allocate(plastic_j2_output(maxval(phase_Noutput), maxNinstance))
plastic_j2_output = ''
allocate(plastic_j2_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(plastic_j2_Noutput(maxNinstance), source=0_pInt)
allocate(plastic_j2_fTaylor(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tau0(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_gdot0(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_n(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_h0(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_h0_slopeLnRate(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_a(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_aTolResistance(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_aTolShear (maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitA(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitB(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitC(maxNinstance), source=0.0_pReal)
allocate(plastic_j2_tausat_SinhFitD(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
phase = phase + 1_pInt ! advance section counter
if (phase_plasticity(phase) == PLASTICITY_J2_ID) then
instance = phase_plasticityInstance(phase)
#ifdef HDF
outID(instance)=HDF5_addGroup(str1,tempResults)
#endif
endif
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_J2_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('flowstress')
plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = flowstress_ID
plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
#ifdef HDF
call HDF5_addScalarDataset(outID(instance),myConstituents,'flowstress','MPa')
allocate(plastic_j2_Output2(instance)%flowstress(myConstituents))
plastic_j2_Output2(instance)%flowstressActive = .true.
#endif
case ('strainrate')
plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = strainrate_ID
plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
#ifdef HDF
call HDF5_addScalarDataset(outID(instance),myConstituents,'strainrate','1/s')
allocate(plastic_j2_Output2(instance)%strainrate(myConstituents))
plastic_j2_Output2(instance)%strainrateActive = .true.
#endif
case default
end select
case ('tau0')
plastic_j2_tau0(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_tau0(instance) < 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('gdot0')
plastic_j2_gdot0(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_gdot0(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('n')
plastic_j2_n(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_n(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('h0')
plastic_j2_h0(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('h0_slope','slopelnrate')
plastic_j2_h0_slopeLnRate(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat')
plastic_j2_tausat(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_tausat(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('tausat_sinhfita')
plastic_j2_tausat_SinhFitA(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitb')
plastic_j2_tausat_SinhFitB(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitc')
plastic_j2_tausat_SinhFitC(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitd')
plastic_j2_tausat_SinhFitD(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('a', 'w0')
plastic_j2_a(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_a(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('taylorfactor')
plastic_j2_fTaylor(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_fTaylor(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('atol_resistance')
plastic_j2_aTolResistance(instance) = IO_floatValue(line,chunkPos,2_pInt)
if (plastic_j2_aTolResistance(instance) <= 0.0_pReal) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
case ('atol_shear')
plastic_j2_aTolShear(instance) = IO_floatValue(line,chunkPos,2_pInt)
case default
end select
endif; endif
enddo parsingFile
initializeInstances: do phase = 1_pInt, size(phase_plasticity)
myPhase: if (phase_plasticity(phase) == PLASTICITY_j2_ID) then
NofMyPhase=count(material_phase==phase)
instance = phase_plasticityInstance(phase)
!--------------------------------------------------------------------------------------------------
! sanity checks
if (plastic_j2_aTolShear(instance) <= 0.0_pReal) &
plastic_j2_aTolShear(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
select case(plastic_j2_outputID(o,instance))
case(flowstress_ID,strainrate_ID)
mySize = 1_pInt
case default
end select
outputFound: if (mySize > 0_pInt) then
plastic_j2_sizePostResult(o,instance) = mySize
plastic_j2_sizePostResults(instance) = &
plastic_j2_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
!--------------------------------------------------------------------------------------------------
! allocate state arrays
sizeState = 2_pInt
sizeDotState = sizeState
sizeDeltaState = 0_pInt
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%sizePostResults = plastic_j2_sizePostResults(instance)
plasticState(phase)%nSlip = 1
plasticState(phase)%nTwin = 0
plasticState(phase)%nTrans= 0
allocate(plasticState(phase)%aTolState ( sizeState))
plasticState(phase)%aTolState(1) = plastic_j2_aTolResistance(instance)
plasticState(phase)%aTolState(2) = plastic_j2_aTolShear(instance)
allocate(plasticState(phase)%state0 ( sizeState,NofMyPhase))
plasticState(phase)%state0(1,1:NofMyPhase) = plastic_j2_tau0(instance)
plasticState(phase)%state0(2,1:NofMyPhase) = 0.0_pReal
allocate(plasticState(phase)%partionedState0 ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%subState0 ( 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)%deltaState (sizeDeltaState,NofMyPhase),source=0.0_pReal)
if (.not. analyticJaco) then
allocate(plasticState(phase)%state_backup ( sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase),source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase),source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase),source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NofMyPhase)
plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NofMyPhase)
endif myPhase
enddo initializeInstances
end subroutine plastic_j2_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_deviatoric33, &
math_mul33xx33
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar_3333 !< derivative of Lp with respect to Tstar as 4th order tensor
real(pReal) :: &
gamma_dot, & !< strainrate
norm_Tstar_dev, & !< euclidean norm of Tstar_dev
squarenorm_Tstar_dev !< square of the euclidean norm of Tstar_dev
integer(pInt) :: &
instance, &
k, l, m, n
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
norm_Tstar_dev = sqrt(squarenorm_Tstar_dev)
if (norm_Tstar_dev <= 0.0_pReal) then ! Tstar == 0 --> both Lp and dLp_dTstar are zero
Lp = 0.0_pReal
dLp_dTstar99 = 0.0_pReal
else
gamma_dot = plastic_j2_gdot0(instance) &
* (sqrt(1.5_pReal) * norm_Tstar_dev / (plastic_j2_fTaylor(instance) * &
plasticState(phaseAt(ipc,ip,el))%state(1,phasememberAt(ipc,ip,el)))) &
**plastic_j2_n(instance)
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/plastic_j2_fTaylor(instance)
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,m,n) = (plastic_j2_n(instance)-1.0_pReal) * &
Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
dLp_dTstar99 = math_Plain3333to99(gamma_dot / plastic_j2_fTaylor(instance) * &
dLp_dTstar_3333 / norm_Tstar_dev)
end if
end subroutine plastic_j2_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_j2_dotState(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
gamma_dot, & !< strainrate
hardening, & !< hardening coefficient
saturation, & !< saturation resistance
norm_Tstar_dev !< euclidean norm of Tstar_dev
integer(pInt) :: &
instance, & !< instance of my instance (unique number of my constitutive model)
of, & !< shortcut notation for offset position in state array
ph !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
of = phasememberAt(ipc,ip,el)
ph = phaseAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!--------------------------------------------------------------------------------------------------
! norm of deviatoric part of 2nd Piola-Kirchhoff stress
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
Tstar_dev_v(4:6) = Tstar_v(4:6)
norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
!--------------------------------------------------------------------------------------------------
! strain rate
gamma_dot = plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!-----------------------------------------------------------------------------------
(plastic_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**plastic_j2_n(instance)
!--------------------------------------------------------------------------------------------------
! hardening coefficient
if (abs(gamma_dot) > 1e-12_pReal) then
if (abs(plastic_j2_tausat_SinhFitA(instance)) <= tiny(0.0_pReal)) then
saturation = plastic_j2_tausat(instance)
else
saturation = ( plastic_j2_tausat(instance) &
+ ( log( ( gamma_dot / plastic_j2_tausat_SinhFitA(instance)&
)**(1.0_pReal / plastic_j2_tausat_SinhFitD(instance))&
+ sqrt( ( gamma_dot / plastic_j2_tausat_SinhFitA(instance) &
)**(2.0_pReal / plastic_j2_tausat_SinhFitD(instance)) &
+ 1.0_pReal ) &
) & ! asinh(K) = ln(K + sqrt(K^2 +1))
)**(1.0_pReal / plastic_j2_tausat_SinhFitC(instance)) &
/ ( plastic_j2_tausat_SinhFitB(instance) &
* (gamma_dot / plastic_j2_gdot0(instance))**(1.0_pReal / plastic_j2_n(instance)) &
) &
)
endif
hardening = ( plastic_j2_h0(instance) + plastic_j2_h0_slopeLnRate(instance) * log(gamma_dot) ) &
* abs( 1.0_pReal - plasticState(ph)%state(1,of)/saturation )**plastic_j2_a(instance) &
* sign(1.0_pReal, 1.0_pReal - plasticState(ph)%state(1,of)/saturation)
else
hardening = 0.0_pReal
endif
plasticState(ph)%dotState(1,of) = hardening * gamma_dot
plasticState(ph)%dotState(2,of) = gamma_dot
end subroutine plastic_j2_dotState
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_j2_postResults(Tstar_v,ipc,ip,el)
use math, only: &
math_mul6x6
use material, only: &
material_phase, &
plasticState, &
phaseAt, phasememberAt, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(plastic_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
plastic_j2_postResults
real(pReal), dimension(6) :: &
Tstar_dev_v ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: &
norm_Tstar_dev ! euclidean norm of Tstar_dev
integer(pInt) :: &
instance, & !< instance of my instance (unique number of my constitutive model)
of, & !< shortcut notation for offset position in state array
ph, & !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
c, &
o
of = phasememberAt(ipc,ip,el)
ph = phaseAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!--------------------------------------------------------------------------------------------------
! calculate deviatoric part of 2nd Piola-Kirchhoff stress and its norm
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
Tstar_dev_v(4:6) = Tstar_v(4:6)
norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
c = 0_pInt
plastic_j2_postResults = 0.0_pReal
outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
select case(plastic_j2_outputID(o,instance))
case (flowstress_ID)
plastic_j2_postResults(c+1_pInt) = plasticState(ph)%state(1,of)
c = c + 1_pInt
case (strainrate_ID)
plastic_j2_postResults(c+1_pInt) = &
plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!----------------------------------------------------------------------------------
(plastic_j2_fTaylor(instance) * plasticState(ph)%state(1,of)) ) ** plastic_j2_n(instance)
c = c + 1_pInt
end select
enddo outputsLoop
end function plastic_j2_postResults
end module plastic_j2

View File

@ -1,109 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for purely elastic material
!--------------------------------------------------------------------------------------------------
module plastic_none
use prec, only: &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
plastic_none_sizePostResults
integer(pInt), dimension(:,:), allocatable, target, public :: &
plastic_none_sizePostResult !< size of each post result output
public :: &
plastic_none_init
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_none_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level, &
debug_constitutive, &
debug_levelBasic
use IO, only: &
IO_timeStamp
use numerics, only: &
worldrank, &
numerics_integrator
use material, only: &
phase_plasticity, &
PLASTICITY_NONE_label, &
material_phase, &
plasticState, &
PLASTICITY_none_ID
implicit none
integer(pInt) :: &
maxNinstance, &
phase, &
NofMyPhase, &
sizeState, &
sizeDotState, &
sizeDeltaState
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONE_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_plasticity == PLASTICITY_none_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
initializeInstances: do phase = 1_pInt, size(phase_plasticity)
if (phase_plasticity(phase) == PLASTICITY_none_ID) then
NofMyPhase=count(material_phase==phase)
sizeState = 0_pInt
plasticState(phase)%sizeState = sizeState
sizeDotState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
sizeDeltaState = 0_pInt
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%sizePostResults = 0_pInt
plasticState(phase)%nSlip = 0_pInt
plasticState(phase)%nTwin = 0_pInt
plasticState(phase)%nTrans = 0_pInt
allocate(plasticState(phase)%aTolState (sizeState))
allocate(plasticState(phase)%state0 (sizeState,NofMyPhase))
allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase))
allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase))
allocate(plasticState(phase)%state (sizeState,NofMyPhase))
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase))
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase))
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase))
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase))
if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase))
allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase))
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase))
if (any(numerics_integrator == 5_pInt)) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase))
endif
enddo initializeInstances
allocate(plastic_none_sizePostResults(maxNinstance), source=0_pInt)
end subroutine plastic_none_init
end module plastic_none

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -1,19 +0,0 @@
# group source
set (SOURCE "source_thermal_dissipation"
"source_thermal_externalheat"
"source_damage_isoBrittle"
"source_damage_isoDuctile"
"source_damage_anisoBrittle"
"source_damage_anisoDuctile"
"source_vacancy_phenoplasticity"
"source_vacancy_irradiation"
"source_vacancy_thermalfluc"
)
# compile module and cumulatively link the
# compiled libraries
foreach (p ${KINEMATICS})
add_library (${p} "${p}.f90")
target_link_libraries(${p} DAMASK_DRIVERS)
add_library (DAMASK_DRIVERS ALIAS ${p})
endforeach (p)

View File

@ -1,425 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Luv Sharma, Max-Planck-Institut fŸr Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut fŸr Eisenforschung GmbH
!> @brief material subroutine incorporating anisotropic brittle damage source mechanism
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_anisoBrittle
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
source_damage_anisoBrittle_sizePostResults, & !< cumulative size of post results
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
source_damage_anisoBrittle_instance !< instance of source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_damage_anisoBrittle_Noutput !< number of outputs per instance of this source
integer(pInt), dimension(:), allocatable, private :: &
source_damage_anisoBrittle_totalNcleavage !< total number of cleavage systems
integer(pInt), dimension(:,:), allocatable, private :: &
source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family
real(pReal), dimension(:), allocatable, private :: &
source_damage_anisoBrittle_aTol, &
source_damage_anisoBrittle_sdot_0, &
source_damage_anisoBrittle_N
real(pReal), dimension(:,:), allocatable, private :: &
source_damage_anisoBrittle_critDisp, &
source_damage_anisoBrittle_critLoad
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
source_damage_anisoBrittle_outputID !< ID of each post result output
public :: &
source_damage_anisoBrittle_init, &
source_damage_anisoBrittle_dotState, &
source_damage_anisobrittle_getRateAndItsTangent, &
source_damage_anisoBrittle_postResults
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_damage_anisoBrittle_label, &
SOURCE_damage_anisoBrittle_ID, &
material_Nphase, &
material_phase, &
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
use lattice, only: &
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoBrittle_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0_pInt)
allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_damage_anisoBrittle_ID) &
source_damage_anisoBrittle_offset(phase) = source
enddo
enddo
allocate(source_damage_anisoBrittle_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance), source=0_pInt)
allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),maxNinstance))
source_damage_anisoBrittle_output = ''
allocate(source_damage_anisoBrittle_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(source_damage_anisoBrittle_Noutput(maxNinstance), source=0_pInt)
allocate(source_damage_anisoBrittle_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(source_damage_anisoBrittle_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt)
allocate(source_damage_anisoBrittle_totalNcleavage(maxNinstance), source=0_pInt)
allocate(source_damage_anisoBrittle_aTol(maxNinstance), source=0.0_pReal)
allocate(source_damage_anisoBrittle_sdot_0(maxNinstance), source=0.0_pReal)
allocate(source_damage_anisoBrittle_N(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_anisoBrittle_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_damage_anisoBrittle_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('anisobrittle_drivingforce')
source_damage_anisoBrittle_Noutput(instance) = source_damage_anisoBrittle_Noutput(instance) + 1_pInt
source_damage_anisoBrittle_outputID(source_damage_anisoBrittle_Noutput(instance),instance) = damage_drivingforce_ID
source_damage_anisoBrittle_output(source_damage_anisoBrittle_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
end select
case ('anisobrittle_atol')
source_damage_anisoBrittle_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('anisobrittle_sdot0')
source_damage_anisoBrittle_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('anisobrittle_ratesensitivity')
source_damage_anisoBrittle_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('ncleavage') !
Nchunks_CleavageFamilies = chunkPos(1) - 1_pInt
do j = 1_pInt, Nchunks_CleavageFamilies
source_damage_anisoBrittle_Ncleavage(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
enddo
case ('anisobrittle_criticaldisplacement')
do j = 1_pInt, Nchunks_CleavageFamilies
source_damage_anisoBrittle_critDisp(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
case ('anisobrittle_criticalload')
do j = 1_pInt, Nchunks_CleavageFamilies
source_damage_anisoBrittle_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
end select
endif; endif
enddo parsingFile
!--------------------------------------------------------------------------------------------------
! sanity checks
sanityChecks: do phase = 1_pInt, material_Nphase
myPhase: if (any(phase_source(:,phase) == SOURCE_damage_anisoBrittle_ID)) then
instance = source_damage_anisoBrittle_instance(phase)
source_damage_anisoBrittle_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,phase),& ! limit active cleavage systems per family to min of available and requested
source_damage_anisoBrittle_Ncleavage(1:lattice_maxNcleavageFamily,instance))
source_damage_anisoBrittle_totalNcleavage(instance) = sum(source_damage_anisoBrittle_Ncleavage(:,instance)) ! how many cleavage systems altogether
if (source_damage_anisoBrittle_aTol(instance) < 0.0_pReal) &
source_damage_anisoBrittle_aTol(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3
if (source_damage_anisoBrittle_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//SOURCE_damage_anisoBrittle_LABEL//')')
if (any(source_damage_anisoBrittle_critDisp(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//SOURCE_damage_anisoBrittle_LABEL//')')
if (any(source_damage_anisoBrittle_critLoad(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//SOURCE_damage_anisoBrittle_LABEL//')')
if (source_damage_anisoBrittle_N(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//SOURCE_damage_anisoBrittle_LABEL//')')
endif myPhase
enddo sanityChecks
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_damage_anisoBrittle_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,source_damage_anisoBrittle_Noutput(instance)
select case(source_damage_anisoBrittle_outputID(o,instance))
case(damage_drivingforce_ID)
mySize = 1_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
source_damage_anisoBrittle_sizePostResult(o,instance) = mySize
source_damage_anisoBrittle_sizePostResults(instance) = source_damage_anisoBrittle_sizePostResults(instance) + mySize
endif
enddo outputsLoop
!--------------------------------------------------------------------------------------------------
! Determine size of state array
sizeDotState = 1_pInt
sizeDeltaState = 0_pInt
sizeState = 1_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoBrittle_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), &
source=source_damage_anisoBrittle_aTol(instance))
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_damage_anisoBrittle_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState, &
material_homog, &
damage, &
damageMapping
use lattice, only: &
lattice_Scleavage_v, &
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
integer(pInt) :: &
phase, &
constituent, &
instance, &
sourceOffset, &
damageOffset, &
homog, &
f, i, index_myFamily
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
homog = material_homog(ip,el)
damageOffset = damageMapping(homog)%p(ip,el)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
do f = 1_pInt,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
traction_crit = source_damage_anisoBrittle_critLoad(f,instance)* &
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
source_damage_anisoBrittle_sdot_0(instance)* &
((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**source_damage_anisoBrittle_N(instance) + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**source_damage_anisoBrittle_N(instance) + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**source_damage_anisoBrittle_N(instance))/ &
source_damage_anisoBrittle_critDisp(f,instance)
enddo
enddo
end subroutine source_damage_anisoBrittle_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer(pInt) :: &
phase, constituent, sourceOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
sourceOffset = source_damage_anisoBrittle_offset(phase)
localphiDot = 1.0_pReal - &
sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
end subroutine source_damage_anisobrittle_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results
!--------------------------------------------------------------------------------------------------
function source_damage_anisoBrittle_postResults(ipc,ip,el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(source_damage_anisoBrittle_sizePostResults( &
source_damage_anisoBrittle_instance(phaseAt(ipc,ip,el)))) :: &
source_damage_anisoBrittle_postResults
integer(pInt) :: &
instance, phase, constituent, sourceOffset, o, c
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
c = 0_pInt
source_damage_anisoBrittle_postResults = 0.0_pReal
do o = 1_pInt,source_damage_anisoBrittle_Noutput(instance)
select case(source_damage_anisoBrittle_outputID(o,instance))
case (damage_drivingforce_ID)
source_damage_anisoBrittle_postResults(c+1_pInt) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1_pInt
end select
enddo
end function source_damage_anisoBrittle_postResults
end module source_damage_anisoBrittle

View File

@ -1,415 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incorporating anisotropic ductile damage source mechanism
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_anisoDuctile
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
source_damage_anisoDuctile_sizePostResults, & !< cumulative size of post results
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_damage_anisoDuctile_Noutput !< number of outputs per instance of this damage
integer(pInt), dimension(:), allocatable, private :: &
source_damage_anisoDuctile_totalNslip !< total number of slip systems
integer(pInt), dimension(:,:), allocatable, private :: &
source_damage_anisoDuctile_Nslip !< number of slip systems per family
real(pReal), dimension(:), allocatable, private :: &
source_damage_anisoDuctile_aTol
real(pReal), dimension(:,:), allocatable, private :: &
source_damage_anisoDuctile_critPlasticStrain
real(pReal), dimension(:), allocatable, private :: &
source_damage_anisoDuctile_sdot_0, &
source_damage_anisoDuctile_N
real(pReal), dimension(:,:), allocatable, private :: &
source_damage_anisoDuctile_critLoad
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
source_damage_anisoDuctile_outputID !< ID of each post result output
public :: &
source_damage_anisoDuctile_init, &
source_damage_anisoDuctile_dotState, &
source_damage_anisoDuctile_getRateAndItsTangent, &
source_damage_anisoDuctile_postResults
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_damage_anisoDuctile_label, &
SOURCE_damage_anisoDuctile_ID, &
material_Nphase, &
material_phase, &
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
use lattice, only: &
lattice_maxNslipFamily, &
lattice_NslipSystem
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_anisoDuctile_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_damage_anisoDuctile_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_damage_anisoDuctile_offset(material_Nphase), source=0_pInt)
allocate(source_damage_anisoDuctile_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_damage_anisoDuctile_ID) &
source_damage_anisoDuctile_offset(phase) = source
enddo
enddo
allocate(source_damage_anisoDuctile_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),maxNinstance))
source_damage_anisoDuctile_output = ''
allocate(source_damage_anisoDuctile_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(source_damage_anisoDuctile_Noutput(maxNinstance), source=0_pInt)
allocate(source_damage_anisoDuctile_critLoad(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal)
allocate(source_damage_anisoDuctile_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal)
allocate(source_damage_anisoDuctile_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt)
allocate(source_damage_anisoDuctile_totalNslip(maxNinstance), source=0_pInt)
allocate(source_damage_anisoDuctile_N(maxNinstance), source=0.0_pReal)
allocate(source_damage_anisoDuctile_sdot_0(maxNinstance), source=0.0_pReal)
allocate(source_damage_anisoDuctile_aTol(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_anisoDuctile_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_damage_anisoDuctile_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('anisoductile_drivingforce')
source_damage_anisoDuctile_Noutput(instance) = source_damage_anisoDuctile_Noutput(instance) + 1_pInt
source_damage_anisoDuctile_outputID(source_damage_anisoDuctile_Noutput(instance),instance) = damage_drivingforce_ID
source_damage_anisoDuctile_output(source_damage_anisoDuctile_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
end select
case ('anisoductile_atol')
source_damage_anisoDuctile_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('nslip') !
Nchunks_SlipFamilies = chunkPos(1) - 1_pInt
do j = 1_pInt, Nchunks_SlipFamilies
source_damage_anisoDuctile_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
enddo
case ('anisoductile_sdot0')
source_damage_anisoDuctile_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('anisoductile_criticalplasticstrain')
do j = 1_pInt, Nchunks_SlipFamilies
source_damage_anisoDuctile_critPlasticStrain(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
case ('anisoductile_ratesensitivity')
source_damage_anisoDuctile_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('anisoductile_criticalload')
do j = 1_pInt, Nchunks_SlipFamilies
source_damage_anisoDuctile_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
enddo
end select
endif; endif
enddo parsingFile
!--------------------------------------------------------------------------------------------------
! sanity checks
sanityChecks: do phase = 1_pInt, size(phase_source)
myPhase: if (any(phase_source(:,phase) == SOURCE_damage_anisoDuctile_ID)) then
instance = source_damage_anisoDuctile_instance(phase)
source_damage_anisoDuctile_Nslip(1:lattice_maxNslipFamily,instance) = &
min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active cleavage systems per family to min of available and requested
source_damage_anisoDuctile_Nslip(1:lattice_maxNslipFamily,instance))
source_damage_anisoDuctile_totalNslip(instance) = sum(source_damage_anisoDuctile_Nslip(:,instance))
if (source_damage_anisoDuctile_aTol(instance) < 0.0_pReal) &
source_damage_anisoDuctile_aTol(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3
if (source_damage_anisoDuctile_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//SOURCE_damage_anisoDuctile_LABEL//')')
if (any(source_damage_anisoDuctile_critPlasticStrain(:,instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//SOURCE_damage_anisoDuctile_LABEL//')')
if (source_damage_anisoDuctile_N(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//SOURCE_damage_anisoDuctile_LABEL//')')
endif myPhase
enddo sanityChecks
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_damage_anisoDuctile_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,source_damage_anisoDuctile_Noutput(instance)
select case(source_damage_anisoDuctile_outputID(o,instance))
case(damage_drivingforce_ID)
mySize = 1_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
source_damage_anisoDuctile_sizePostResult(o,instance) = mySize
source_damage_anisoDuctile_sizePostResults(instance) = source_damage_anisoDuctile_sizePostResults(instance) + mySize
endif
enddo outputsLoop
!--------------------------------------------------------------------------------------------------
! Determine size of state array
sizeDotState = 1_pInt
sizeDeltaState = 0_pInt
sizeState = 1_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoDuctile_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), &
source=source_damage_anisoDuctile_aTol(instance))
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_damage_anisoDuctile_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
sourceState, &
material_homog, &
damage, &
damageMapping
use lattice, only: &
lattice_maxNslipFamily
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer(pInt) :: &
phase, &
constituent, &
sourceOffset, &
homog, damageOffset, &
instance, &
index, f, i
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
homog = material_homog(ip,el)
damageOffset = damageMapping(homog)%p(ip,el)
index = 1_pInt
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
do f = 1_pInt,lattice_maxNslipFamily
do i = 1_pInt,source_damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
plasticState(phase)%slipRate(index,constituent)/ &
((damage(homog)%p(damageOffset))**source_damage_anisoDuctile_N(instance))/ &
source_damage_anisoDuctile_critPlasticStrain(f,instance)
index = index + 1_pInt
enddo
enddo
end subroutine source_damage_anisoDuctile_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer(pInt) :: &
phase, constituent, sourceOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
sourceOffset = source_damage_anisoDuctile_offset(phase)
localphiDot = 1.0_pReal - &
sourceState(phase)%p(sourceOffset)%state(1,constituent)* &
phi
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
end subroutine source_damage_anisoDuctile_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results
!--------------------------------------------------------------------------------------------------
function source_damage_anisoDuctile_postResults(ipc,ip,el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(source_damage_anisoDuctile_sizePostResults( &
source_damage_anisoDuctile_instance(phaseAt(ipc,ip,el)))) :: &
source_damage_anisoDuctile_postResults
integer(pInt) :: &
instance, phase, constituent, sourceOffset, o, c
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
c = 0_pInt
source_damage_anisoDuctile_postResults = 0.0_pReal
do o = 1_pInt,source_damage_anisoDuctile_Noutput(instance)
select case(source_damage_anisoDuctile_outputID(o,instance))
case (damage_drivingforce_ID)
source_damage_anisoDuctile_postResults(c+1_pInt) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1_pInt
end select
enddo
end function source_damage_anisoDuctile_postResults
end module source_damage_anisoDuctile

View File

@ -1,383 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incoprorating isotropic brittle damage source mechanism
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_isoBrittle
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
source_damage_isoBrittle_sizePostResults, & !< cumulative size of post results
source_damage_isoBrittle_offset, & !< which source is my current damage mechanism?
source_damage_isoBrittle_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
source_damage_isoBrittle_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_isoBrittle_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_damage_isoBrittle_Noutput !< number of outputs per instance of this damage
real(pReal), dimension(:), allocatable, private :: &
source_damage_isoBrittle_aTol, &
source_damage_isoBrittle_N, &
source_damage_isoBrittle_critStrainEnergy
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 ToDo
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
source_damage_isoBrittle_outputID !< ID of each post result output
public :: &
source_damage_isoBrittle_init, &
source_damage_isoBrittle_deltaState, &
source_damage_isoBrittle_getRateAndItsTangent, &
source_damage_isoBrittle_postResults
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_damage_isoBrittle_label, &
SOURCE_damage_isoBrittle_ID, &
material_Nphase, &
material_phase, &
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoBrittle_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_damage_isoBrittle_offset(material_Nphase), source=0_pInt)
allocate(source_damage_isoBrittle_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
source_damage_isoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoBrittle_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_damage_isoBrittle_ID) &
source_damage_isoBrittle_offset(phase) = source
enddo
enddo
allocate(source_damage_isoBrittle_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),maxNinstance))
source_damage_isoBrittle_output = ''
allocate(source_damage_isoBrittle_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(source_damage_isoBrittle_Noutput(maxNinstance), source=0_pInt)
allocate(source_damage_isoBrittle_critStrainEnergy(maxNinstance), source=0.0_pReal)
allocate(source_damage_isoBrittle_N(maxNinstance), source=1.0_pReal)
allocate(source_damage_isoBrittle_aTol(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_isoBrittle_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_damage_isoBrittle_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('isobrittle_drivingforce')
source_damage_isoBrittle_Noutput(instance) = source_damage_isoBrittle_Noutput(instance) + 1_pInt
source_damage_isoBrittle_outputID(source_damage_isoBrittle_Noutput(instance),instance) = damage_drivingforce_ID
source_damage_isoBrittle_output(source_damage_isoBrittle_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
end select
case ('isobrittle_criticalstrainenergy')
source_damage_isoBrittle_critStrainEnergy(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('isobrittle_n')
source_damage_isoBrittle_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('isobrittle_atol')
source_damage_isoBrittle_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif; endif
enddo parsingFile
!--------------------------------------------------------------------------------------------------
! sanity checks
sanityChecks: do phase = 1_pInt, material_Nphase
myPhase: if (any(phase_source(:,phase) == SOURCE_damage_isoBrittle_ID)) then
instance = source_damage_isoBrittle_instance(phase)
if (source_damage_isoBrittle_aTol(instance) < 0.0_pReal) &
source_damage_isoBrittle_aTol(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3
if (source_damage_isoBrittle_critStrainEnergy(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='criticalStrainEnergy ('//SOURCE_damage_isoBrittle_LABEL//')')
endif myPhase
enddo sanityChecks
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_damage_isoBrittle_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase)
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,source_damage_isoBrittle_Noutput(instance)
select case(source_damage_isoBrittle_outputID(o,instance))
case(damage_drivingforce_ID)
mySize = 1_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
source_damage_isoBrittle_sizePostResult(o,instance) = mySize
source_damage_isoBrittle_sizePostResults(instance) = source_damage_isoBrittle_sizePostResults(instance) + mySize
endif
enddo outputsLoop
! Determine size of state array
sizeDotState = 1_pInt
sizeDeltaState = 1_pInt
sizeState = 1_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoBrittle_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), &
source=source_damage_isoBrittle_aTol(instance))
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_damage_isoBrittle_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState, &
material_homog, &
phase_NstiffnessDegradations, &
phase_stiffnessDegradation, &
porosity, &
porosityMapping, &
STIFFNESS_DEGRADATION_porosity_ID
use math, only : &
math_mul33x33, &
math_mul66x6, &
math_Mandel33to6, &
math_transpose33, &
math_I3
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
Fe
real(pReal), intent(in), dimension(6,6) :: &
C
integer(pInt) :: &
phase, constituent, instance, sourceOffset, mech
real(pReal) :: &
strain(6), &
stiffness(6,6), &
strainenergy
phase = phaseAt(ipc,ip,el) !< phase ID at ipc,ip,el
constituent = phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el
! ToDo: capability for multiple instances of SAME source within given phase. Needs Ninstance loop from here on!
instance = source_damage_isoBrittle_instance(phase) !< instance of damage_isoBrittle source
sourceOffset = source_damage_isoBrittle_offset(phase)
stiffness = C
do mech = 1_pInt, phase_NstiffnessDegradations(phase)
select case(phase_stiffnessDegradation(mech,phase))
case (STIFFNESS_DEGRADATION_porosity_ID)
stiffness = porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* &
porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* &
stiffness
end select
enddo
strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
strainenergy = 2.0_pReal*sum(strain*math_mul66x6(stiffness,strain))/ &
source_damage_isoBrittle_critStrainEnergy(instance)
if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent)
else
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
endif
end subroutine source_damage_isoBrittle_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer(pInt) :: &
phase, constituent, instance, sourceOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase)
localphiDot = (1.0_pReal - phi)**(source_damage_isoBrittle_N(instance) - 1.0_pReal) - &
phi*sourceState(phase)%p(sourceOffset)%state(1,constituent)
dLocalphiDot_dPhi = - (source_damage_isoBrittle_N(instance) - 1.0_pReal)* &
(1.0_pReal - phi)**max(0.0_pReal,source_damage_isoBrittle_N(instance) - 2.0_pReal) &
- sourceState(phase)%p(sourceOffset)%state(1,constituent)
end subroutine source_damage_isoBrittle_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results
!--------------------------------------------------------------------------------------------------
function source_damage_isoBrittle_postResults(ipc,ip,el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(source_damage_isoBrittle_sizePostResults( &
source_damage_isoBrittle_instance(phaseAt(ipc,ip,el)))) :: &
source_damage_isoBrittle_postResults
integer(pInt) :: &
instance, phase, constituent, sourceOffset, o, c
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase)
c = 0_pInt
source_damage_isoBrittle_postResults = 0.0_pReal
do o = 1_pInt,source_damage_isoBrittle_Noutput(instance)
select case(source_damage_isoBrittle_outputID(o,instance))
case (damage_drivingforce_ID)
source_damage_isoBrittle_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
end select
enddo
end function source_damage_isoBrittle_postResults
end module source_damage_isoBrittle

View File

@ -1,350 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut fŸr Eisenforschung GmbH
!> @author Luv Sharma, Max-Planck-Institut fŸr Eisenforschung GmbH
!> @brief material subroutine incoprorating isotropic ductile damage source mechanism
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_isoDuctile
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
source_damage_isoDuctile_sizePostResults, & !< cumulative size of post results
source_damage_isoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_isoDuctile_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
source_damage_isoDuctile_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_isoDuctile_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_damage_isoDuctile_Noutput !< number of outputs per instance of this damage
real(pReal), dimension(:), allocatable, private :: &
source_damage_isoDuctile_aTol, &
source_damage_isoDuctile_critPlasticStrain, &
source_damage_isoDuctile_N
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 ToDo
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
source_damage_isoDuctile_outputID !< ID of each post result output
public :: &
source_damage_isoDuctile_init, &
source_damage_isoDuctile_dotState, &
source_damage_isoDuctile_getRateAndItsTangent, &
source_damage_isoDuctile_postResults
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_damage_isoDuctile_label, &
SOURCE_damage_isoDuctile_ID, &
material_Nphase, &
material_phase, &
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_damage_isoDuctile_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_damage_isoDuctile_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_damage_isoDuctile_offset(material_Nphase), source=0_pInt)
allocate(source_damage_isoDuctile_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
source_damage_isoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoDuctile_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_damage_isoDuctile_ID) &
source_damage_isoDuctile_offset(phase) = source
enddo
enddo
allocate(source_damage_isoDuctile_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),maxNinstance))
source_damage_isoDuctile_output = ''
allocate(source_damage_isoDuctile_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(source_damage_isoDuctile_Noutput(maxNinstance), source=0_pInt)
allocate(source_damage_isoDuctile_critPlasticStrain(maxNinstance), source=0.0_pReal)
allocate(source_damage_isoDuctile_N(maxNinstance), source=0.0_pReal)
allocate(source_damage_isoDuctile_aTol(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_isoDuctile_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_damage_isoDuctile_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('isoductile_drivingforce')
source_damage_isoDuctile_Noutput(instance) = source_damage_isoDuctile_Noutput(instance) + 1_pInt
source_damage_isoDuctile_outputID(source_damage_isoDuctile_Noutput(instance),instance) = damage_drivingforce_ID
source_damage_isoDuctile_output(source_damage_isoDuctile_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
end select
case ('isoductile_criticalplasticstrain')
source_damage_isoDuctile_critPlasticStrain(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('isoductile_ratesensitivity')
source_damage_isoDuctile_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('isoductile_atol')
source_damage_isoDuctile_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif; endif
enddo parsingFile
!--------------------------------------------------------------------------------------------------
! sanity checks
sanityChecks: do phase = 1_pInt, material_Nphase
myPhase: if (any(phase_source(:,phase) == SOURCE_damage_isoDuctile_ID)) then
instance = source_damage_isoDuctile_instance(phase)
if (source_damage_isoDuctile_aTol(instance) < 0.0_pReal) &
source_damage_isoDuctile_aTol(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3
if (source_damage_isoDuctile_critPlasticStrain(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='critical plastic strain ('//SOURCE_damage_isoDuctile_LABEL//')')
endif myPhase
enddo sanityChecks
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_damage_isoDuctile_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,source_damage_isoDuctile_Noutput(instance)
select case(source_damage_isoDuctile_outputID(o,instance))
case(damage_drivingforce_ID)
mySize = 1_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
source_damage_isoDuctile_sizePostResult(o,instance) = mySize
source_damage_isoDuctile_sizePostResults(instance) = source_damage_isoDuctile_sizePostResults(instance) + mySize
endif
enddo outputsLoop
! Determine size of state array
sizeDotState = 1_pInt
sizeDeltaState = 0_pInt
sizeState = 1_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoDuctile_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), &
source=source_damage_isoDuctile_aTol(instance))
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_damage_isoDuctile_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
sourceState, &
material_homog, &
damage, &
damageMapping
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer(pInt) :: &
phase, constituent, instance, homog, sourceOffset, damageOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
homog = material_homog(ip,el)
damageOffset = damageMapping(homog)%p(ip,el)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sum(plasticState(phase)%slipRate(:,constituent))/ &
((damage(homog)%p(damageOffset))**source_damage_isoDuctile_N(instance))/ &
source_damage_isoDuctile_critPlasticStrain(instance)
end subroutine source_damage_isoDuctile_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer(pInt) :: &
phase, constituent, sourceOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
sourceOffset = source_damage_isoDuctile_offset(phase)
localphiDot = 1.0_pReal - &
sourceState(phase)%p(sourceOffset)%state(1,constituent)* &
phi
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
end subroutine source_damage_isoDuctile_getRateAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of local damage results
!--------------------------------------------------------------------------------------------------
function source_damage_isoDuctile_postResults(ipc,ip,el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), dimension(source_damage_isoDuctile_sizePostResults( &
source_damage_isoDuctile_instance(phaseAt(ipc,ip,el)))) :: &
source_damage_isoDuctile_postResults
integer(pInt) :: &
instance, phase, constituent, sourceOffset, o, c
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
c = 0_pInt
source_damage_isoDuctile_postResults = 0.0_pReal
do o = 1_pInt,source_damage_isoDuctile_Noutput(instance)
select case(source_damage_isoDuctile_outputID(o,instance))
case (damage_drivingforce_ID)
source_damage_isoDuctile_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
end select
enddo
end function source_damage_isoDuctile_postResults
end module source_damage_isoDuctile

View File

@ -1,220 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for thermal source due to plastic dissipation
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_thermal_dissipation
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
source_thermal_dissipation_sizePostResults, & !< cumulative size of post results
source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
source_thermal_dissipation_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_thermal_dissipation_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_thermal_dissipation_Noutput !< number of outputs per instance of this source
real(pReal), dimension(:), allocatable, private :: &
source_thermal_dissipation_coldworkCoeff
public :: &
source_thermal_dissipation_init, &
source_thermal_dissipation_getRateAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_dissipation_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_thermal_dissipation_label, &
SOURCE_thermal_dissipation_ID, &
material_Nphase, &
material_phase, &
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_thermal_dissipation_offset(material_Nphase), source=0_pInt)
allocate(source_thermal_dissipation_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
source_thermal_dissipation_instance(phase) = count(phase_source(:,1:phase) == SOURCE_thermal_dissipation_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == SOURCE_thermal_dissipation_ID) &
source_thermal_dissipation_offset(phase) = source
enddo
enddo
allocate(source_thermal_dissipation_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_thermal_dissipation_output (maxval(phase_Noutput),maxNinstance))
source_thermal_dissipation_output = ''
allocate(source_thermal_dissipation_Noutput(maxNinstance), source=0_pInt)
allocate(source_thermal_dissipation_coldworkCoeff(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_thermal_dissipation_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_thermal_dissipation_instance(phase) ! which instance of my source is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('dissipation_coldworkcoeff')
source_thermal_dissipation_coldworkCoeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif; endif
enddo parsingFile
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_thermal_dissipation_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_thermal_dissipation_instance(phase)
sourceOffset = source_thermal_dissipation_offset(phase)
sizeDotState = 0_pInt
sizeDeltaState = 0_pInt
sizeState = 0_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_thermal_dissipation_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_thermal_dissipation_init
!--------------------------------------------------------------------------------------------------
!> @brief returns local vacancy generation rate
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar_v, Lp, ipc, ip, el)
use math, only: &
math_Mandel6to33
use material, only: &
phaseAt, phasememberAt
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal), intent(in), dimension(3,3) :: &
Lp
real(pReal), intent(out) :: &
TDot, &
dTDOT_dT
integer(pInt) :: &
instance, phase, constituent
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_thermal_dissipation_instance(phase)
TDot = source_thermal_dissipation_coldworkCoeff(instance)* &
sum(abs(math_Mandel6to33(Tstar_v)*Lp))
dTDOT_dT = 0.0_pReal
end subroutine source_thermal_dissipation_getRateAndItsTangent
end module source_thermal_dissipation

View File

@ -1,277 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for thermal source due to plastic dissipation
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_thermal_externalheat
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
source_thermal_externalheat_sizePostResults, & !< cumulative size of post results
source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
source_thermal_externalheat_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_thermal_externalheat_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_thermal_externalheat_Noutput !< number of outputs per instance of this source
integer(pInt), dimension(:), allocatable, private :: &
source_thermal_externalheat_nIntervals
real(pReal), dimension(:,:), allocatable, private :: &
source_thermal_externalheat_time, &
source_thermal_externalheat_rate
public :: &
source_thermal_externalheat_init, &
source_thermal_externalheat_dotState, &
source_thermal_externalheat_getRateAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_thermal_externalheat_label, &
SOURCE_thermal_externalheat_ID, &
material_Nphase, &
material_phase, &
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase,interval
character(len=65536) :: &
tag = '', &
line = ''
real(pReal), allocatable, dimension(:,:) :: temp_time, temp_rate
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_thermal_externalheat_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_thermal_externalheat_offset(material_Nphase), source=0_pInt)
allocate(source_thermal_externalheat_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
source_thermal_externalheat_instance(phase) = count(phase_source(:,1:phase) == SOURCE_thermal_externalheat_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == SOURCE_thermal_externalheat_ID) &
source_thermal_externalheat_offset(phase) = source
enddo
enddo
allocate(source_thermal_externalheat_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_thermal_externalheat_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_thermal_externalheat_output (maxval(phase_Noutput),maxNinstance))
source_thermal_externalheat_output = ''
allocate(source_thermal_externalheat_Noutput(maxNinstance), source=0_pInt)
allocate(source_thermal_externalheat_nIntervals(maxNinstance), source=0_pInt)
allocate(temp_time(maxNinstance,1000), source=0.0_pReal)
allocate(temp_rate(maxNinstance,1000), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_thermal_externalheat_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_thermal_externalheat_instance(phase) ! which instance of my source is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('externalheat_time')
if (chunkPos(1) <= 2_pInt) &
call IO_error(150_pInt,ext_msg=trim(tag)//' ('//SOURCE_thermal_externalheat_label//')')
source_thermal_externalheat_nIntervals(instance) = chunkPos(1) - 2_pInt
do interval = 1, source_thermal_externalheat_nIntervals(instance) + 1_pInt
temp_time(instance, interval) = IO_floatValue(line,chunkPos,1_pInt + interval)
enddo
case ('externalheat_rate')
do interval = 1, source_thermal_externalheat_nIntervals(instance) + 1_pInt
temp_rate(instance, interval) = IO_floatValue(line,chunkPos,1_pInt + interval)
enddo
end select
endif; endif
enddo parsingFile
allocate(source_thermal_externalheat_time(maxNinstance,maxval(source_thermal_externalheat_nIntervals)+1_pInt), source=0.0_pReal)
allocate(source_thermal_externalheat_rate(maxNinstance,maxval(source_thermal_externalheat_nIntervals)+1_pInt), source=0.0_pReal)
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_thermal_externalheat_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_thermal_externalheat_instance(phase)
sourceOffset = source_thermal_externalheat_offset(phase)
source_thermal_externalheat_time(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt) = &
temp_time(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt)
source_thermal_externalheat_rate(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt) = &
temp_rate(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt)
sizeDotState = 1_pInt
sizeDeltaState = 0_pInt
sizeState = 1_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_thermal_externalheat_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.00001_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_thermal_externalheat_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_dotState(ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer(pInt) :: &
phase, &
constituent, &
sourceOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
sourceOffset = source_thermal_externalheat_offset(phase)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 1.0_pReal
end subroutine source_thermal_externalheat_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local vacancy generation rate
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out) :: &
TDot, &
dTDot_dT
integer(pInt) :: &
instance, phase, constituent, sourceOffset, interval
real(pReal) :: &
norm_time
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_thermal_externalheat_instance(phase)
sourceOffset = source_thermal_externalheat_offset(phase)
do interval = 1, source_thermal_externalheat_nIntervals(instance)
norm_time = (sourceState(phase)%p(sourceOffset)%state(1,constituent) - &
source_thermal_externalheat_time(instance,interval)) / &
(source_thermal_externalheat_time(instance,interval+1) - &
source_thermal_externalheat_time(instance,interval))
if (norm_time >= 0.0_pReal .and. norm_time < 1.0_pReal) &
TDot = source_thermal_externalheat_rate(instance,interval ) * (1.0_pReal - norm_time) + &
source_thermal_externalheat_rate(instance,interval+1) * norm_time
enddo
dTDot_dT = 0.0
end subroutine source_thermal_externalheat_getRateAndItsTangent
end module source_thermal_externalheat

View File

@ -1,253 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for vacancy generation due to irradiation
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_vacancy_irradiation
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
source_vacancy_irradiation_sizePostResults, & !< cumulative size of post results
source_vacancy_irradiation_offset, & !< which source is my current damage mechanism?
source_vacancy_irradiation_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
source_vacancy_irradiation_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_vacancy_irradiation_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_vacancy_irradiation_Noutput !< number of outputs per instance of this damage
real(pReal), dimension(:), allocatable, private :: &
source_vacancy_irradiation_cascadeProb, &
source_vacancy_irradiation_cascadeVolume
public :: &
source_vacancy_irradiation_init, &
source_vacancy_irradiation_deltaState, &
source_vacancy_irradiation_getRateAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_vacancy_irradiation_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_vacancy_irradiation_label, &
SOURCE_vacancy_irradiation_ID, &
material_Nphase, &
material_phase, &
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_irradiation_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_vacancy_irradiation_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_vacancy_irradiation_offset(material_Nphase), source=0_pInt)
allocate(source_vacancy_irradiation_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
source_vacancy_irradiation_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_irradiation_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_vacancy_irradiation_ID) &
source_vacancy_irradiation_offset(phase) = source
enddo
enddo
allocate(source_vacancy_irradiation_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_vacancy_irradiation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_vacancy_irradiation_output(maxval(phase_Noutput),maxNinstance))
source_vacancy_irradiation_output = ''
allocate(source_vacancy_irradiation_Noutput(maxNinstance), source=0_pInt)
allocate(source_vacancy_irradiation_cascadeProb(maxNinstance), source=0.0_pReal)
allocate(source_vacancy_irradiation_cascadeVolume(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_irradiation_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_vacancy_irradiation_instance(phase) ! which instance of my vacancy is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('irradiation_cascadeprobability')
source_vacancy_irradiation_cascadeProb(instance) = IO_floatValue(line,chunkPos,2_pInt)
case ('irradiation_cascadevolume')
source_vacancy_irradiation_cascadeVolume(instance) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif; endif
enddo parsingFile
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_vacancy_irradiation_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_vacancy_irradiation_instance(phase)
sourceOffset = source_vacancy_irradiation_offset(phase)
sizeDotState = 2_pInt
sizeDeltaState = 2_pInt
sizeState = 2_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_irradiation_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.1_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_vacancy_irradiation_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_vacancy_irradiation_deltaState(ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer(pInt) :: &
phase, constituent, sourceOffset
real(pReal) :: &
randNo
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
sourceOffset = source_vacancy_irradiation_offset(phase)
call random_number(randNo)
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
randNo - sourceState(phase)%p(sourceOffset)%state(1,constituent)
call random_number(randNo)
sourceState(phase)%p(sourceOffset)%deltaState(2,constituent) = &
randNo - sourceState(phase)%p(sourceOffset)%state(2,constituent)
end subroutine source_vacancy_irradiation_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns local vacancy generation rate
!--------------------------------------------------------------------------------------------------
subroutine source_vacancy_irradiation_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out) :: &
CvDot, dCvDot_dCv
integer(pInt) :: &
instance, phase, constituent, sourceOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_vacancy_irradiation_instance(phase)
sourceOffset = source_vacancy_irradiation_offset(phase)
CvDot = 0.0_pReal
dCvDot_dCv = 0.0_pReal
if (sourceState(phase)%p(sourceOffset)%state0(1,constituent) < source_vacancy_irradiation_cascadeProb(instance)) &
CvDot = sourceState(phase)%p(sourceOffset)%state0(2,constituent)*source_vacancy_irradiation_cascadeVolume(instance)
end subroutine source_vacancy_irradiation_getRateAndItsTangent
end module source_vacancy_irradiation

View File

@ -1,215 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for vacancy generation due to plasticity
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_vacancy_phenoplasticity
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
source_vacancy_phenoplasticity_sizePostResults, & !< cumulative size of post results
source_vacancy_phenoplasticity_offset, & !< which source is my current damage mechanism?
source_vacancy_phenoplasticity_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
source_vacancy_phenoplasticity_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_vacancy_phenoplasticity_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_vacancy_phenoplasticity_Noutput !< number of outputs per instance of this damage
real(pReal), dimension(:), allocatable, private :: &
source_vacancy_phenoplasticity_rateCoeff
public :: &
source_vacancy_phenoplasticity_init, &
source_vacancy_phenoplasticity_getRateAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_vacancy_phenoplasticity_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_vacancy_phenoplasticity_label, &
SOURCE_vacancy_phenoplasticity_ID, &
material_Nphase, &
material_phase, &
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_phenoplasticity_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_vacancy_phenoplasticity_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_vacancy_phenoplasticity_offset(material_Nphase), source=0_pInt)
allocate(source_vacancy_phenoplasticity_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
source_vacancy_phenoplasticity_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_phenoplasticity_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_vacancy_phenoplasticity_ID) &
source_vacancy_phenoplasticity_offset(phase) = source
enddo
enddo
allocate(source_vacancy_phenoplasticity_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_vacancy_phenoplasticity_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_vacancy_phenoplasticity_output(maxval(phase_Noutput),maxNinstance))
source_vacancy_phenoplasticity_output = ''
allocate(source_vacancy_phenoplasticity_Noutput(maxNinstance), source=0_pInt)
allocate(source_vacancy_phenoplasticity_rateCoeff(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_phenoplasticity_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_vacancy_phenoplasticity_instance(phase) ! which instance of my vacancy is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('phenoplasticity_ratecoeff')
source_vacancy_phenoplasticity_rateCoeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif; endif
enddo parsingFile
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_vacancy_phenoplasticity_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_vacancy_phenoplasticity_instance(phase)
sourceOffset = source_vacancy_phenoplasticity_offset(phase)
sizeDotState = 0_pInt
sizeDeltaState = 0_pInt
sizeState = 0_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_phenoplasticity_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_vacancy_phenoplasticity_init
!--------------------------------------------------------------------------------------------------
!> @brief returns local vacancy generation rate
!--------------------------------------------------------------------------------------------------
subroutine source_vacancy_phenoplasticity_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
plasticState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out) :: &
CvDot, dCvDot_dCv
integer(pInt) :: &
instance, phase, constituent
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_vacancy_phenoplasticity_instance(phase)
CvDot = &
source_vacancy_phenoplasticity_rateCoeff(instance)* &
sum(plasticState(phase)%slipRate(:,constituent))
dCvDot_dCv = 0.0_pReal
end subroutine source_vacancy_phenoplasticity_getRateAndItsTangent
end module source_vacancy_phenoplasticity

View File

@ -1,255 +0,0 @@
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for vacancy generation due to thermal fluctuations
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_vacancy_thermalfluc
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
source_vacancy_thermalfluc_sizePostResults, & !< cumulative size of post results
source_vacancy_thermalfluc_offset, & !< which source is my current damage mechanism?
source_vacancy_thermalfluc_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
source_vacancy_thermalfluc_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_vacancy_thermalfluc_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_vacancy_thermalfluc_Noutput !< number of outputs per instance of this damage
real(pReal), dimension(:), allocatable, private :: &
source_vacancy_thermalfluc_amplitude, &
source_vacancy_thermalfluc_normVacancyEnergy
public :: &
source_vacancy_thermalfluc_init, &
source_vacancy_thermalfluc_deltaState, &
source_vacancy_thermalfluc_getRateAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_vacancy_thermalfluc_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use lattice, only: &
lattice_vacancyFormationEnergy
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_vacancy_thermalfluc_label, &
SOURCE_vacancy_thermalfluc_ID, &
material_Nphase, &
material_phase, &
sourceState, &
MATERIAL_partPhase
use numerics,only: &
analyticJaco, &
worldrank, &
numerics_integrator
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
character(len=65536) :: &
tag = '', &
line = ''
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_thermalfluc_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_source == SOURCE_vacancy_thermalfluc_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_vacancy_thermalfluc_offset(material_Nphase), source=0_pInt)
allocate(source_vacancy_thermalfluc_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
source_vacancy_thermalfluc_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_thermalfluc_ID)
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == source_vacancy_thermalfluc_ID) &
source_vacancy_thermalfluc_offset(phase) = source
enddo
enddo
allocate(source_vacancy_thermalfluc_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_vacancy_thermalfluc_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_vacancy_thermalfluc_output(maxval(phase_Noutput),maxNinstance))
source_vacancy_thermalfluc_output = ''
allocate(source_vacancy_thermalfluc_Noutput(maxNinstance), source=0_pInt)
allocate(source_vacancy_thermalfluc_amplitude(maxNinstance), source=0.0_pReal)
allocate(source_vacancy_thermalfluc_normVacancyEnergy(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_thermalfluc_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_vacancy_thermalfluc_instance(phase) ! which instance of my vacancy is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('thermalfluctuation_amplitude')
source_vacancy_thermalfluc_amplitude(instance) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif; endif
enddo parsingFile
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_vacancy_thermalfluc_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_vacancy_thermalfluc_instance(phase)
source_vacancy_thermalfluc_normVacancyEnergy(instance) = &
lattice_vacancyFormationEnergy(phase)/1.3806488e-23_pReal
sourceOffset = source_vacancy_thermalfluc_offset(phase)
sizeDotState = 1_pInt
sizeDeltaState = 1_pInt
sizeState = 1_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_thermalfluc_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.1_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (.not. analyticJaco) then
allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_vacancy_thermalfluc_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_vacancy_thermalfluc_deltaState(ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer(pInt) :: &
phase, constituent, sourceOffset
real(pReal) :: &
randNo
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
sourceOffset = source_vacancy_thermalfluc_offset(phase)
call random_number(randNo)
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
randNo - 0.5_pReal - sourceState(phase)%p(sourceOffset)%state(1,constituent)
end subroutine source_vacancy_thermalfluc_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns local vacancy generation rate
!--------------------------------------------------------------------------------------------------
subroutine source_vacancy_thermalfluc_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
material_homog, &
temperature, &
thermalMapping, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out) :: &
CvDot, dCvDot_dCv
integer(pInt) :: &
instance, phase, constituent, sourceOffset
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_vacancy_thermalfluc_instance(phase)
sourceOffset = source_vacancy_thermalfluc_offset(phase)
CvDot = source_vacancy_thermalfluc_amplitude(instance)* &
sourceState(phase)%p(sourceOffset)%state0(2,constituent)* &
exp(-source_vacancy_thermalfluc_normVacancyEnergy(instance)/ &
temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el)))
dCvDot_dCv = 0.0_pReal
end subroutine source_vacancy_thermalfluc_getRateAndItsTangent
end module source_vacancy_thermalfluc