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

Merge remote-tracking branch 'remotes/origin/clean-and-polish-damage' into development

This commit is contained in:
Franz Roters 2020-03-13 07:19:52 +01:00
parent e93141de3a d6e9eb71ea
commit e8c16e7f98
28 changed files with 1349 additions and 1691 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 9573ce7bd2c1a7188c1aac5b83aa76d480c2bdb0
Subproject commit d0d5b5a22be9778187b100214c782747793bb956

2
examples/ConfigFiles/Phase_Isotropic_AluminumIsotropic.config
View File

@ -7,7 +7,7 @@ plasticity isotropic
(output) flowstress
(output) strainrate
lattice_structure isotropic
lattice_structure iso
c11 110.9e9
c12 58.34e9

2
examples/ConfigFiles/Phase_Isotropic_FreeSurface.config
View File

@ -12,7 +12,7 @@ plasticity isotropic
(output) flowstress
(output) strainrate
lattice_structure isotropic
lattice_structure iso
c11 10e9
c12 0.0
gdot0 0.001

3
examples/ConfigFiles/Phase_None_IsotropicVolumePreservation.config
View File

@ -1,9 +1,8 @@
[IsotropicVolumePreservation]
elasticity hooke
plasticity none
### Material parameters ###
lattice_structure isotropic
lattice_structure iso
C11 100.0e9
C12 66.6666667e9

4
examples/SpectralMethod/EshelbyInclusion/material.config
View File

@ -16,7 +16,7 @@ t0 330.0
#.................
[isotropic matrix]
lattice_structure isotropic
lattice_structure iso
plasticity none
{config/elastic_isotropic.config}
{config/thermal.config}
@ -45,7 +45,7 @@ plasticity none
#.................
[isotropic inclusion]
lattice_structure isotropic
lattice_structure iso
plasticity none
{config/elastic_isotropic.config}
{config/thermal.config}

21
src/CPFEM.f90
View File

@ -27,9 +27,6 @@ module CPFEM
implicit none
private
real(pReal), parameter, private :: &
CPFEM_odd_stress = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle
CPFEM_odd_jacobian = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle
real(pReal), dimension (:,:,:), allocatable, private :: &
CPFEM_cs !< Cauchy stress
real(pReal), dimension (:,:,:,:), allocatable, private :: &
@ -150,7 +147,10 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog, mySource
logical updateJaco ! flag indicating if JAcobian has to be updated
logical updateJaco ! flag indicating if Jacobian has to be updated
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle
elCP = mesh_FEM2DAMASK_elem(elFE)
@ -193,8 +193,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress
CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6)
CPFEM_cs(1:6,ip,elCP) = rnd * ODD_STRESS
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
chosenThermal2: select case (thermal_type(material_homogenizationAt(elCP)))
case (THERMAL_conduction_ID) chosenThermal2
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
@ -226,8 +226,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
endif
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = rnd*CPFEM_odd_stress
CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian*math_identity2nd(6)
CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
!*** deformation gradient is not outdated
@ -257,8 +257,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress
CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6)
CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
else terminalIllness
@ -331,6 +331,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
end subroutine CPFEM_general
!--------------------------------------------------------------------------------------------------
!> @brief Forward data for new time increment.
!--------------------------------------------------------------------------------------------------

2
src/DAMASK_marc.f90
View File

@ -40,7 +40,7 @@ module DAMASK_interface
implicit none
private
logical, public :: symmetricSolver
logical, protected, public :: symmetricSolver
character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat'
public :: &

14
src/HDF5_utilities.f90
View File

@ -77,20 +77,6 @@ module HDF5_utilities
module procedure HDF5_addAttribute_real_array
end interface HDF5_addAttribute
!--------------------------------------------------------------------------------------------------
public :: &
HDF5_utilities_init, &
HDF5_openFile, &
HDF5_closeFile, &
HDF5_addAttribute, &
HDF5_closeGroup ,&
HDF5_openGroup, &
HDF5_addGroup, &
HDF5_read, &
HDF5_write, &
HDF5_setLink, &
HDF5_objectExists
contains

8
src/constitutive_plastic_dislotwin.f90
View File

@ -198,10 +198,10 @@ module subroutine plastic_dislotwin_init
prm%n0_sl = lattice_slip_normal(prm%N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID) &
prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == lattice_FCC_ID) &
.and. (prm%N_sl(1) == 12)
if(prm%fccTwinTransNucleation) &
prm%fcc_twinNucleationSlipPair = lattice_fcc_twinNucleationSlipPair
prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR
prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl))
prm%rho_dip_0 = config%getFloats('rhoedgedip0',requiredSize=size(prm%N_sl))
@ -230,7 +230,7 @@ module subroutine plastic_dislotwin_init
prm%omega = config%getFloat('omega', defaultVal = 1000.0_pReal) &
* merge(12.0_pReal, &
8.0_pReal, &
lattice_structure(p) == LATTICE_FCC_ID .or. lattice_structure(p) == LATTICE_HEX_ID)
lattice_structure(p) == lattice_FCC_ID .or. lattice_structure(p) == lattice_HEX_ID)
! expand: family => system
@ -335,7 +335,7 @@ module subroutine plastic_dislotwin_init
config%getFloat('a_bcc', defaultVal=0.0_pReal), &
config%getFloat('a_fcc', defaultVal=0.0_pReal))
if (lattice_structure(p) /= LATTICE_fcc_ID) then
if (lattice_structure(p) /= lattice_FCC_ID) then
prm%dot_N_0_tr = config%getFloats('ndot0_trans')
prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr)
endif

42
src/damage_local.f90
View File

@ -16,15 +16,9 @@ module damage_local
implicit none
private
enum, bind(c)
enumerator :: &
undefined_ID, &
damage_ID
end enum
type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tparameters), dimension(:), allocatable :: &
@ -33,7 +27,7 @@ module damage_local
public :: &
damage_local_init, &
damage_local_updateState, &
damage_local_Results
damage_local_results
contains
@ -43,29 +37,18 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine damage_local_init
integer :: maxNinstance,o,NofMyHomog,h
character(len=pStringLen), dimension(:), allocatable :: outputs
integer :: Ninstance,NofMyHomog,h
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'; flush(6)
maxNinstance = count(damage_type == DAMAGE_local_ID)
if (maxNinstance == 0) return
Ninstance = count(damage_type == DAMAGE_local_ID)
allocate(param(Ninstance))
allocate(param(maxNinstance))
do h = 1, size(damage_type)
do h = 1, size(config_homogenization)
if (damage_type(h) /= DAMAGE_LOCAL_ID) cycle
associate(prm => param(damage_typeInstance(h)),config => config_homogenization(h))
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do o=1, size(outputs)
select case(outputs(o))
case ('damage')
prm%outputID = [prm%outputID , damage_ID]
end select
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
NofMyHomog = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1
@ -179,14 +162,13 @@ subroutine damage_local_results(homog,group)
integer, intent(in) :: homog
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(damage_typeInstance(homog)))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
outputsLoop: do o = 1,size(prm%output)
select case(prm%output(o))
case ('damage')
call results_writeDataset(group,damage(homog)%p,'phi',&
'damage indicator','-')
end select

60
src/damage_nonlocal.f90
View File

@ -18,15 +18,9 @@ module damage_nonlocal
implicit none
private
enum, bind(c)
enumerator :: &
undefined_ID, &
damage_ID
end enum
type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tparameters), dimension(:), allocatable :: &
@ -35,10 +29,10 @@ module damage_nonlocal
public :: &
damage_nonlocal_init, &
damage_nonlocal_getSourceAndItsTangent, &
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getDiffusion, &
damage_nonlocal_getMobility, &
damage_nonlocal_putNonLocalDamage, &
damage_nonlocal_Results
damage_nonlocal_results
contains
@ -48,29 +42,18 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_init
integer :: maxNinstance,o,NofMyHomog,h
character(len=pStringLen), dimension(:), allocatable :: outputs
integer :: Ninstance,NofMyHomog,h
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'; flush(6)
maxNinstance = count(damage_type == DAMAGE_nonlocal_ID)
if (maxNinstance == 0) return
Ninstance = count(damage_type == DAMAGE_nonlocal_ID)
allocate(param(Ninstance))
allocate(param(maxNinstance))
do h = 1, size(damage_type)
do h = 1, size(config_homogenization)
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
associate(prm => param(damage_typeInstance(h)),config => config_homogenization(h))
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do o=1, size(outputs)
select case(outputs(o))
case ('damage')
prm%outputID = [prm%outputID, damage_ID]
end select
enddo
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
NofMyHomog = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1
@ -145,28 +128,28 @@ end subroutine damage_nonlocal_getSourceAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized non local damage diffusion tensor in reference configuration
!--------------------------------------------------------------------------------------------------
function damage_nonlocal_getDiffusion33(ip,el)
function damage_nonlocal_getDiffusion(ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
damage_nonlocal_getDiffusion33
damage_nonlocal_getDiffusion
integer :: &
homog, &
grain
homog = material_homogenizationAt(el)
damage_nonlocal_getDiffusion33 = 0.0_pReal
damage_nonlocal_getDiffusion = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog)
damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + &
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phaseAt(grain,el)))
damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion(1:3,1:3,material_phaseAt(grain,el)))
enddo
damage_nonlocal_getDiffusion33 = &
charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
damage_nonlocal_getDiffusion = &
charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Ngrains(homog),pReal)
end function damage_nonlocal_getDiffusion33
end function damage_nonlocal_getDiffusion
!--------------------------------------------------------------------------------------------------
@ -220,14 +203,13 @@ subroutine damage_nonlocal_results(homog,group)
integer, intent(in) :: homog
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(damage_typeInstance(homog)))
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
outputsLoop: do o = 1,size(prm%output)
select case(prm%output(o))
case ('damage')
call results_writeDataset(group,damage(homog)%p,'phi',&
'damage indicator','-')
end select

4
src/grid/grid_damage_spectral.f90
View File

@ -249,7 +249,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion33(1,cell) - K_ref, &
vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(1,cell) - K_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward
@ -294,7 +294,7 @@ subroutine updateReference
mu_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
K_ref = K_ref + damage_nonlocal_getDiffusion33(1,cell)
K_ref = K_ref + damage_nonlocal_getDiffusion(1,cell)
mu_ref = mu_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo
K_ref = K_ref*wgt

7
src/grid/grid_thermal_spectral.f90
View File

@ -252,7 +252,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity33(1,cell) - K_ref, &
vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(1,cell) - K_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward
@ -294,9 +294,8 @@ subroutine updateReference
mu_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
K_ref = K_ref + thermal_conduction_getConductivity33(1,cell)
mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)
K_ref = K_ref + thermal_conduction_getConductivity(1,cell)
mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* thermal_conduction_getSpecificHeat(1,cell)
enddo; enddo; enddo
K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)

166
src/kinematics_cleavage_opening.f90
View File

@ -27,25 +27,12 @@ module kinematics_cleavage_opening
sdot0, &
n
real(pReal), dimension(:), allocatable :: &
critDisp, &
critLoad
end type
real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems
end type tParameters
! Begin Deprecated
integer, dimension(:), allocatable :: &
kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems
integer, dimension(:,:), allocatable :: &
kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family
real(pReal), dimension(:), allocatable :: &
kinematics_cleavage_opening_sdot_0, &
kinematics_cleavage_opening_N
real(pReal), dimension(:,:), allocatable :: &
kinematics_cleavage_opening_critDisp, &
kinematics_cleavage_opening_critLoad
! End Deprecated
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
kinematics_cleavage_opening_init, &
@ -60,61 +47,54 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_init
integer, allocatable, dimension(:) :: tempInt
real(pReal), allocatable, dimension(:) :: tempFloat
integer :: maxNinstance,p,instance
integer :: Ninstance,p
character(len=pStringLen) :: extmsg = ''
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'; flush(6)
maxNinstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID)
if (maxNinstance == 0) return
Ninstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0)
do p = 1, size(config_phase)
kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct?
enddo
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)
allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0)
allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal)
allocate(param(Ninstance))
do p = 1, size(config_phase)
kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID)
if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle
instance = kinematics_cleavage_opening_instance(p)
kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0')
kinematics_cleavage_opening_N(instance) = config_phase(p)%getFloat('anisobrittle_ratesensitivity')
tempInt = config_phase(p)%getInts('ncleavage')
kinematics_cleavage_opening_Ncleavage(1:size(tempInt),instance) = tempInt
tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt))
kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat
associate(prm => param(kinematics_cleavage_opening_instance(p)), &
config => config_phase(p))
tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt))
kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat
prm%Ncleavage = config%getInts('ncleavage')
prm%totalNcleavage = sum(prm%Ncleavage)
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,p),& ! 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,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
prm%n = config%getFloat('anisobrittle_ratesensitivity')
prm%sdot0 = config%getFloat('anisobrittle_sdot0')
prm%critLoad = config%getFloats('anisobrittle_criticalload',requiredSize=size(prm%Ncleavage))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(prm%Ncleavage,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system
prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage)
! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_n'
if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critLoad'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
end associate
enddo
end subroutine kinematics_cleavage_opening_init
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
@ -130,72 +110,54 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
integer :: &
instance, phase, &
homog, damageOffset, &
f, i, index_myFamily, k, l, m, n
i, k, l, m, n
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
phase = material_phaseAt(ipc,el)
instance = kinematics_cleavage_opening_instance(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal
do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
traction_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)
associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el))))
do i = 1,prm%totalNcleavage
traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset)**2.0_pReal
traction_d = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,1,i))
if (abs(traction_d) > traction_crit + tol_math_check) then
udotd = sign(1.0_pReal,traction_d)* prm%sdot0 * ((abs(traction_d) - traction_crit)/traction_crit)**prm%n
Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i)
dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%n / (abs(traction_d) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,1,index_myFamily+i,phase)
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i)
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)
traction_t = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,2,i))
if (abs(traction_t) > traction_crit + tol_math_check) then
udott = sign(1.0_pReal,traction_t)* prm%sdot0 * ((abs(traction_t) - traction_crit)/traction_crit)**prm%n
Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i)
dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%n / (abs(traction_t) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,2,index_myFamily+i,phase)
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i)
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)
traction_n = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,3,i))
if (abs(traction_n) > traction_crit + tol_math_check) then
udotn = sign(1.0_pReal,traction_n)* prm%sdot0 * ((abs(traction_n) - traction_crit)/traction_crit)**prm%n
Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i)
dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%n / (abs(traction_n) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,3,index_myFamily+i,phase)
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i)
endif
enddo
enddo
end associate
end subroutine kinematics_cleavage_opening_LiAndItsTangent

107
src/kinematics_slipplane_opening.f90
View File

@ -28,10 +28,10 @@ module kinematics_slipplane_opening
n
real(pReal), dimension(:), allocatable :: &
critLoad
real(pReal), dimension(:,:), allocatable :: &
slip_direction, &
slip_normal, &
slip_transverse
real(pReal), dimension(:,:,:), allocatable :: &
P_d, &
P_t, &
P_n
end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
@ -49,56 +49,64 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_init
integer :: maxNinstance,p,instance
integer :: Ninstance,p,i
character(len=pStringLen) :: extmsg = ''
real(pReal), dimension(:,:), allocatable :: d,n,t
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'; flush(6)
maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID)
if (maxNinstance == 0) return
Ninstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0)
do p = 1, size(config_phase)
kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct?
enddo
allocate(param(maxNinstance))
allocate(param(Ninstance))
do p = 1, size(config_phase)
kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID)
if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle
associate(prm => param(kinematics_slipplane_opening_instance(p)), &
config => config_phase(p))
instance = kinematics_slipplane_opening_instance(p)
prm%sdot0 = config_phase(p)%getFloat('anisoductile_sdot0')
prm%n = config_phase(p)%getFloat('anisoductile_ratesensitivity')
prm%sdot0 = config%getFloat('anisoductile_sdot0')
prm%n = config%getFloat('anisoductile_ratesensitivity')
prm%Nslip = config%getInts('nslip')
prm%totalNslip = sum(prm%Nslip)
prm%critLoad = config_phase(p)%getFloats('anisoductile_criticalload',requiredSize=size(prm%Nslip ))
d = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
t = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
n = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
allocate(prm%P_d(3,3,size(d,2)),prm%P_t(3,3,size(t,2)),prm%P_n(3,3,size(n,2)))
do i=1, size(n,2)
prm%P_d(1:3,1:3,i) = math_outer(d(1:3,i), n(1:3,i))
prm%P_t(1:3,1:3,i) = math_outer(t(1:3,i), n(1:3,i))
prm%P_n(1:3,1:3,i) = math_outer(n(1:3,i), n(1:3,i))
enddo
prm%critLoad = config%getFloats('anisoductile_criticalload',requiredSize=size(prm%Nslip))
! expand: family => system
prm%critLoad = math_expand(prm%critLoad, prm%Nslip)
prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_n'
if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_sdot0'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad'
! if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) &
! call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) &
! call IO_error(211,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) &
! call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')')
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')')
end associate
enddo
end subroutine kinematics_slipplane_opening_init
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
@ -114,8 +122,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< 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
integer :: &
instance, phase, &
homog, damageOffset, &
@ -134,13 +141,9 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
dLd_dTstar = 0.0_pReal
do i = 1, prm%totalNslip
projection_d = math_outer(prm%slip_direction(1:3,i), prm%slip_normal(1:3,i))
projection_t = math_outer(prm%slip_transverse(1:3,i),prm%slip_normal(1:3,i))
projection_n = math_outer(prm%slip_normal(1:3,i), prm%slip_normal(1:3,i))
traction_d = math_mul33xx33(S,projection_d)
traction_t = math_mul33xx33(S,projection_t)
traction_n = math_mul33xx33(S,projection_n)
traction_d = math_tensordot(S,prm%P_d(1:3,1:3,i))
traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i))
traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i))
traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage
@ -151,18 +154,32 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
udotn = prm%sdot0* ( max(0.0_pReal,traction_n)/traction_crit &
- max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n
if (dNeq0(traction_d)) then
dudotd_dt = udotd*prm%n/traction_d
else
dudotd_dt = 0.0_pReal
endif
if (dNeq0(traction_t)) then
dudott_dt = udott*prm%n/traction_t
else
dudott_dt = 0.0_pReal
endif
if (dNeq0(traction_n)) then
dudotn_dt = udotn*prm%n/traction_n
else
dudotn_dt = 0.0_pReal
endif
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotd_dt*projection_d(k,l)*projection_d(m,n) &
+ dudott_dt*projection_t(k,l)*projection_t(m,n) &
+ dudotn_dt*projection_n(k,l)*projection_n(m,n)
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudotd_dt*prm%P_d(k,l,i)*prm%P_d(m,n,i) &
+ dudott_dt*prm%P_t(k,l,i)*prm%P_t(m,n,i) &
+ dudotn_dt*prm%P_n(k,l,i)*prm%P_n(m,n,i)
Ld = Ld + udotd*projection_d &
+ udott*projection_t &
+ udotn*projection_n
Ld = Ld &
+ udotd*prm%P_d(1:3,1:3,i) &
+ udott*prm%P_t(1:3,1:3,i) &
+ udotn*prm%P_n(1:3,1:3,i)
enddo
end associate

85
src/kinematics_thermal_expansion.f90
View File

@ -15,9 +15,13 @@ module kinematics_thermal_expansion
implicit none
private
integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance
type :: tParameters
real(pReal), allocatable, dimension(:,:,:) :: &
expansion
real(pReal) :: &
T_ref
real(pReal), dimension(3,3,3) :: &
expansion = 0.0_pReal
end type tParameters
type(tParameters), dimension(:), allocatable :: param
@ -36,35 +40,39 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_init
integer :: &
Ninstance, &
p, i
real(pReal), dimension(:), allocatable :: &
temp
integer :: Ninstance,p,i
real(pReal), dimension(:), allocatable :: temp
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'; flush(6)
Ninstance = count(phase_kinematics == KINEMATICS_thermal_expansion_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(kinematics_thermal_expansion_instance(size(config_phase)), source=0)
allocate(param(Ninstance))
do p = 1, size(phase_kinematics)
do p = 1, size(config_phase)
kinematics_thermal_expansion_instance(p) = count(phase_kinematics(:,1:p) == KINEMATICS_thermal_expansion_ID)
if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle
! ToDo: Here we need to decide how to extend the concept of instances to
! kinetics and sources. I would suggest that the same mechanism exists at maximum once per phase
associate(prm => param(kinematics_thermal_expansion_instance(p)), &
config => config_phase(p))
prm%T_ref = config%getFloat('reference_temperature', defaultVal=0.0_pReal)
! read up to three parameters (constant, linear, quadratic with T)
temp = config_phase(p)%getFloats('thermal_expansion11')
!lattice_thermalExpansion33(1,1,1:size(temp),p) = temp
temp = config_phase(p)%getFloats('thermal_expansion22', &
defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
!lattice_thermalExpansion33(2,2,1:size(temp),p) = temp
temp = config_phase(p)%getFloats('thermal_expansion33', &
defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
temp = config%getFloats('thermal_expansion11')
prm%expansion(1,1,1:size(temp)) = temp
temp = config%getFloats('thermal_expansion22',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
prm%expansion(2,2,1:size(temp)) = temp
temp = config%getFloats('thermal_expansion33',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
prm%expansion(3,3,1:size(temp)) = temp
do i=1, size(prm%expansion,3)
prm%expansion(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%expansion(1:3,1:3,i),config%getString('lattice_structure'))
enddo
end associate
enddo
end subroutine kinematics_thermal_expansion_init
@ -77,18 +85,18 @@ pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset)
integer, intent(in) :: &
phase, &
homog, offset
homog, &
offset
real(pReal), dimension(3,3) :: &
kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though)
associate(prm => param(kinematics_thermal_expansion_instance(phase)))
kinematics_thermal_expansion_initialStrain = &
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * &
lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * &
lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * &
lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient
(temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient
(temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient
(temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient
end associate
end function kinematics_thermal_expansion_initialStrain
@ -106,29 +114,30 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip,
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
integer :: &
phase, &
homog, offset
homog
real(pReal) :: &
T, TRef, TDot
T, TDot
phase = material_phaseAt(ipc,el)
homog = material_homogenizationAt(el)
offset = thermalMapping(homog)%p(ip,el)
T = temperature(homog)%p(offset)
TDot = temperatureRate(homog)%p(offset)
TRef = lattice_referenceTemperature(phase)
T = temperature(homog)%p(thermalMapping(homog)%p(ip,el))
TDot = temperatureRate(homog)%p(thermalMapping(homog)%p(ip,el))
associate(prm => param(kinematics_thermal_expansion_instance(phase)))
Li = TDot * ( &
lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient
prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
) / &
(1.0_pReal &
+ lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. &
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. &
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. &
+ prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+ prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+ prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
)
end associate
dLi_dTstar = 0.0_pReal
end subroutine kinematics_thermal_expansion_LiAndItsTangent

1034
src/lattice.f90
View File

File diff suppressed because it is too large Load Diff

33
src/prec.f90
View File

@ -11,30 +11,31 @@ module prec
implicit none
public
! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds
integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit)
integer, parameter :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit)
#if(INT==8)
integer, parameter, public :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit)
integer, parameter :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit)
#else
integer, parameter, public :: pInt = selected_int_kind(9) !< number with at least up to +-1e9 (typically 32 bit)
integer, parameter :: pInt = selected_int_kind(9) !< number with at least up to +-1e9 (typically 32 bit)
#endif
integer, parameter, public :: pLongInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit)
integer, parameter, public :: pStringLen = 256 !< default string length
integer, parameter, public :: pPathLen = 4096 !< maximum length of a path name on linux
integer, parameter :: pLongInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit)
integer, parameter :: pStringLen = 256 !< default string length
integer, parameter :: pPathLen = 4096 !< maximum length of a path name on linux
real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation)
real(pReal), parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation)
type, public :: group_float !< variable length datatype used for storage of state
type :: group_float !< variable length datatype used for storage of state
real(pReal), dimension(:), pointer :: p
end type group_float
type, public :: group_int
type :: group_int
integer, dimension(:), pointer :: p
end type group_int
! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array
type, public :: tState
type :: tState
integer :: &
sizeState = 0, & !< size of state
sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
@ -57,7 +58,7 @@ module prec
RKCK45dotState
end type
type, extends(tState), public :: tPlasticState
type, extends(tState) :: tPlasticState
logical :: &
nonlocal = .false.
real(pReal), pointer, dimension(:,:) :: &
@ -65,22 +66,22 @@ module prec
accumulatedSlip !< accumulated plastic slip
end type
type, public :: tSourceState
type :: tSourceState
type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase
end type
type, public :: tHomogMapping
type :: tHomogMapping
integer, pointer, dimension(:,:) :: p
end type
real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal) !< minimum positive number such that 1.0 + EPSILON /= 1.0.
real(pReal), private, parameter :: PREAL_MIN = tiny(0.0_pReal) !< smallest normalized floating point number
integer, dimension(0), parameter, public :: &
integer, dimension(0), parameter :: &
emptyIntArray = [integer::]
real(pReal), dimension(0), parameter, public :: &
real(pReal), dimension(0), parameter :: &
emptyRealArray = [real(pReal)::]
character(len=pStringLen), dimension(0), parameter, public :: &
character(len=pStringLen), dimension(0), parameter :: &
emptyStringArray = [character(len=pStringLen)::]
private :: &

11
src/quaternions.f90
View File

@ -10,7 +10,7 @@ module quaternions
use IO
implicit none
public
private
real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion.
@ -96,8 +96,12 @@ module quaternions
module procedure aimag__
end interface aimag
private :: &
unitTest
public :: &
quaternions_init, &
assignment(=), &
conjg, aimag, &
log, exp, &
real
contains
@ -511,6 +515,7 @@ subroutine unitTest
q_2 = conjg(q_2) - inverse(q_2)
if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(0,ext_msg='inverse/conjg')
endif
if(dNeq(dot_product(qu,qu),dot_product(q,q))) call IO_error(0,ext_msg='dot_product')
#if !(defined(__GFORTRAN__) && __GNUC__ < 9)
if (norm2(aimag(q)) > 0.0_pReal) then

160
src/source_damage_anisoBrittle.f90
View File

@ -9,8 +9,8 @@ module source_damage_anisoBrittle
use debug
use IO
use math
use material
use discretization
use material
use config
use lattice
use results
@ -22,15 +22,6 @@ module source_damage_anisoBrittle
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
source_damage_anisoBrittle_instance !< instance of source mechanism
integer, dimension(:,:), allocatable :: &
source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum
type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
aTol, &
@ -45,8 +36,8 @@ module source_damage_anisoBrittle
totalNcleavage
integer, dimension(:), allocatable :: &
Ncleavage
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
@ -67,56 +58,39 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_init
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p ,i
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
integer :: Ninstance,sourceOffset,NofMyPhase,p
character(len=pStringLen) :: extmsg = ''
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'; flush(6)
Ninstance = count(phase_source == SOURCE_damage_anisoBrittle_ID)
if (Ninstance == 0) return
Ninstance = count(phase_source == SOURCE_DAMAGE_ANISOBRITTLE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0)
allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0)
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_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0)
allocate(source_damage_anisoBrittle_offset (size(config_phase)), source=0)
allocate(source_damage_anisoBrittle_instance(size(config_phase)), source=0)
allocate(param(Ninstance))
do p = 1, size(config_phase)
source_damage_anisoBrittle_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ANISOBRITTLE_ID)
do sourceOffset = 1, phase_Nsources(p)
if (phase_source(sourceOffset,p) == SOURCE_DAMAGE_ANISOBRITTLE_ID) then
source_damage_anisoBrittle_offset(p) = sourceOffset
exit
endif
enddo
if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISOBRITTLE_ID)) cycle
associate(prm => param(source_damage_anisoBrittle_instance(p)), &
config => config_phase(p))
prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal)
prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray)
prm%totalNcleavage = sum(prm%Ncleavage)
prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal)
prm%N = config%getFloat('anisobrittle_ratesensitivity')
prm%sdot_0 = config%getFloat('anisobrittle_sdot0')
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity'
if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0'
prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray)
prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage))
prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage))
@ -127,41 +101,26 @@ subroutine source_damage_anisoBrittle_init
prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage)
prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage)
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload'
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement'
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_n'
if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critLoad'
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critDisp'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
case ('anisobrittle_drivingforce')
prm%outputID = [prm%outputID, damage_drivingforce_ID]
end select
enddo
NofMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
end associate
phase = p
NofMyPhase=count(material_phaseAt==phase) * discretization_nIP
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
source_damage_anisoBrittle_Ncleavage(1:size(param(instance)%Ncleavage),instance) = param(instance)%Ncleavage
enddo
end subroutine source_damage_anisoBrittle_init
@ -178,49 +137,42 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S
integer :: &
phase, &
constituent, &
instance, &
sourceOffset, &
damageOffset, &
homog, &
f, i, index_myFamily, index
i
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit
phase = material_phaseAt(ipc,el)
constituent = material_phasememberAt(ipc,ip,el)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
associate(prm => param(source_damage_anisoBrittle_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
do i = 1, prm%totalNcleavage
index = 1
do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,1,i))
traction_t = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,2,i))
traction_n = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,3,i))
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
traction_crit = prm%critLoad(i)*damage(homog)%p(damageOffset)**2.0_pReal
traction_crit = param(instance)%critLoad(index)* &
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
+ prm%sdot_0 / prm%critDisp(i) &
* ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%N + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%N + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%N)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
param(instance)%sdot_0* &
((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**param(instance)%N + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**param(instance)%N + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ &
param(instance)%critDisp(index)
index = index + 1
enddo
enddo
end associate
end subroutine source_damage_anisoBrittle_dotState
@ -238,16 +190,17 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
sourceOffset
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)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
end subroutine source_damage_anisoBrittle_getRateAndItsTangent
@ -258,15 +211,14 @@ subroutine source_damage_anisoBrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
integer :: sourceOffset, o, instance
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
integer :: o
associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_drivingforce_ID)
associate(prm => param(source_damage_anisoBrittle_instance(phase)), &
stt => sourceState(phase)%p(source_damage_anisoBrittle_offset(phase))%state)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('anisobrittle_drivingforce')
call results_writeDataset(group,stt,'tbd','driving force','tbd')
end select
enddo outputsLoop

121
src/source_damage_anisoDuctile.f90
View File

@ -21,12 +21,6 @@ module source_damage_anisoDuctile
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
aTol, &
@ -37,8 +31,8 @@ module source_damage_anisoDuctile
totalNslip
integer, dimension(:), allocatable :: &
Nslip
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
@ -59,58 +53,46 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_init
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p ,i
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
integer :: Ninstance,sourceOffset,NofMyPhase,p
character(len=pStringLen) :: extmsg = ''
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'; flush(6)
Ninstance = count(phase_source == SOURCE_damage_anisoDuctile_ID)
if (Ninstance == 0) return
Ninstance = count(phase_source == SOURCE_DAMAGE_ANISODUCTILE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoDuctile_offset (size(config_phase)), source=0)
allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0)
do phase = 1, size(config_phase)
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(param(Ninstance))
do p = 1, size(config_phase)
source_damage_anisoDuctile_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ANISODUCTILE_ID)
do sourceOffset = 1, phase_Nsources(p)
if (phase_source(sourceOffset,p) == SOURCE_DAMAGE_ANISODUCTILE_ID) then
source_damage_anisoDuctile_offset(p) = sourceOffset
exit
endif
enddo
if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISODUCTILE_ID)) cycle
associate(prm => param(source_damage_anisoDuctile_instance(p)), &
config => config_phase(p))
prm%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal)
prm%N = config%getFloat('anisoductile_ratesensitivity')
prm%totalNslip = sum(prm%Nslip)
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity'
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
prm%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal)
prm%N = config%getFloat('anisoductile_ratesensitivity')
prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip))
! expand: family => system
prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip)
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity'
if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain'
!--------------------------------------------------------------------------------------------------
@ -119,30 +101,13 @@ subroutine source_damage_anisoDuctile_init
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
case ('anisoductile_drivingforce')
prm%outputID = [prm%outputID, damage_drivingforce_ID]
end select
enddo
NofMyPhase=count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
end associate
phase = p
NofMyPhase=count(material_phaseAt==phase) * discretization_nIP
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
enddo
end subroutine source_damage_anisoDuctile_init
@ -157,28 +122,28 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: &
phase, &
constituent, &
sourceOffset, &
homog, damageOffset, &
instance, &
damageOffset, &
homog, &
i
phase = material_phaseAt(ipc,el)
constituent = material_phasememberAt(ipc,ip,el)
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
do i = 1, param(instance)%totalNslip
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
plasticState(phase)%slipRate(i,constituent)/ &
((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i)
associate(prm => param(source_damage_anisoDuctile_instance(phase)))
do i = 1, prm%totalNslip
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
+ plasticState(phase)%slipRate(i,constituent)/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain(i)
enddo
end associate
end subroutine source_damage_anisoDuctile_dotState
@ -196,16 +161,17 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
sourceOffset
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)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
end subroutine source_damage_anisoDuctile_getRateAndItsTangent
@ -216,15 +182,14 @@ subroutine source_damage_anisoDuctile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
integer :: sourceOffset, o, instance
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
integer :: o
associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_drivingforce_ID)
associate(prm => param(source_damage_anisoDuctile_instance(phase)), &
stt => sourceState(phase)%p(source_damage_anisoDuctile_offset(phase))%state)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('anisoductile_drivingforce')
call results_writeDataset(group,stt,'tbd','driving force','tbd')
end select
enddo outputsLoop

121
src/source_damage_isoBrittle.f90
View File

@ -16,24 +16,18 @@ module source_damage_isoBrittle
implicit none
private
integer, dimension(:), allocatable :: &
source_damage_isoBrittle_offset, &
source_damage_isoBrittle_instance
enum, bind(c)
enumerator :: &
undefined_ID, &
damage_drivingforce_ID
end enum
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
critStrainEnergy, &
N, &
aTol
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
@ -43,7 +37,7 @@ module source_damage_isoBrittle
source_damage_isoBrittle_init, &
source_damage_isoBrittle_deltaState, &
source_damage_isoBrittle_getRateAndItsTangent, &
source_damage_isoBrittle_Results
source_damage_isoBrittle_results
contains
@ -54,49 +48,38 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_init
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p,i
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
integer :: Ninstance,sourceOffset,NofMyPhase,p
character(len=pStringLen) :: extmsg = ''
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>'; flush(6)
Ninstance = count(phase_source == SOURCE_damage_isoBrittle_ID)
if (Ninstance == 0) return
Ninstance = count(phase_source == SOURCE_DAMAGE_ISOBRITTLE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_isoBrittle_offset(material_Nphase), source=0)
allocate(source_damage_isoBrittle_instance(material_Nphase), source=0)
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_offset (size(config_phase)), source=0)
allocate(source_damage_isoBrittle_instance(size(config_phase)), source=0)
allocate(param(Ninstance))
do p = 1, size(config_phase)
source_damage_isoBrittle_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ISOBRITTLE_ID)
do sourceOffset = 1, phase_Nsources(p)
if (phase_source(sourceOffset,p) == SOURCE_DAMAGE_ISOBRITTLE_ID) then
source_damage_isoBrittle_offset(p) = sourceOffset
exit
endif
enddo
if (all(phase_source(:,p) /= SOURCE_DAMAGE_ISOBRITTLE_ID)) cycle
associate(prm => param(source_damage_isoBrittle_instance(p)), &
config => config_phase(p))
prm%aTol = config%getFloat('isobrittle_atol',defaultVal = 1.0e-3_pReal)
prm%N = config%getFloat('isobrittle_n')
prm%critStrainEnergy = config%getFloat('isobrittle_criticalstrainenergy')
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_atol'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_n'
if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy'
@ -107,34 +90,18 @@ subroutine source_damage_isoBrittle_init
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
case ('isobrittle_drivingforce')
prm%outputID = [prm%outputID, damage_drivingforce_ID]
end select
enddo
NofMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,1)
sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
end associate
phase = p
NofMyPhase = count(material_phaseAt==phase) * discretization_nIP
instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,1)
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
enddo
end subroutine source_damage_isoBrittle_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
@ -148,22 +115,24 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
Fe
real(pReal), intent(in), dimension(6,6) :: &
C
integer :: &
phase, constituent, instance, sourceOffset
phase, &
constituent, &
sourceOffset
real(pReal), dimension(6) :: &
strain
real(pReal) :: &
strain(6), &
strainenergy
phase = material_phaseAt(ipc,el) !< phase ID at ipc,ip,el
constituent = material_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)
strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/param(instance)%critStrainEnergy
associate(prm => param(source_damage_isoBrittle_instance(phase)))
strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%critStrainEnergy
! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/param(instance)%critStrainEnergy
if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
@ -174,9 +143,11 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
endif
end associate
end subroutine source_damage_isoBrittle_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
@ -190,17 +161,18 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
instance, sourceOffset
instance = source_damage_isoBrittle_instance(phase)
integer :: &
sourceOffset
sourceOffset = source_damage_isoBrittle_offset(phase)
localphiDot = (1.0_pReal - phi)**(param(instance)%N - 1.0_pReal) - &
phi*sourceState(phase)%p(sourceOffset)%state(1,constituent)
dLocalphiDot_dPhi = - (param(instance)%N - 1.0_pReal)* &
(1.0_pReal - phi)**max(0.0_pReal,param(instance)%N - 2.0_pReal) &
associate(prm => param(source_damage_isoBrittle_instance(phase)))
localphiDot = (1.0_pReal - phi)**(prm%n - 1.0_pReal) &
- phi*sourceState(phase)%p(sourceOffset)%state(1,constituent)
dLocalphiDot_dPhi = - (prm%n - 1.0_pReal)* (1.0_pReal - phi)**max(0.0_pReal,prm%n - 2.0_pReal) &
- sourceState(phase)%p(sourceOffset)%state(1,constituent)
end associate
end subroutine source_damage_isoBrittle_getRateAndItsTangent
@ -212,15 +184,14 @@ subroutine source_damage_isoBrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
integer :: sourceOffset, o, instance
instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase)
integer :: o
associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_drivingforce_ID)
associate(prm => param(source_damage_isoBrittle_instance(phase)), &
stt => sourceState(phase)%p(source_damage_isoBrittle_offset(phase))%state)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('isobrittle_drivingforce')
call results_writeDataset(group,stt,'tbd','driving force','tbd')
end select
enddo outputsLoop

110
src/source_damage_isoDuctile.f90
View File

@ -15,22 +15,18 @@ module source_damage_isoDuctile
implicit none
private
integer, dimension(:), allocatable :: &
source_damage_isoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_isoDuctile_instance !< instance of damage source mechanism
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ToDo
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
critPlasticStrain, &
N, &
aTol
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
@ -51,49 +47,38 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_init
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p,i
integer(kind(undefined_ID)) :: &
outputID
integer :: Ninstance,sourceOffset,NofMyPhase,p
character(len=pStringLen) :: extmsg = ''
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init -+>>>'
Ninstance = count(phase_source == SOURCE_damage_isoDuctile_ID)
if (Ninstance == 0) return
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init -+>>>'; flush(6)
Ninstance = count(phase_source == SOURCE_DAMAGE_ISODUCTILE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_isoDuctile_offset(material_Nphase), source=0)
allocate(source_damage_isoDuctile_instance(material_Nphase), source=0)
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_offset (size(config_phase)), source=0)
allocate(source_damage_isoDuctile_instance(size(config_phase)), source=0)
allocate(param(Ninstance))
do p = 1, size(config_phase)
source_damage_isoDuctile_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ISODUCTILE_ID)
do sourceOffset = 1, phase_Nsources(p)
if (phase_source(sourceOffset,p) == SOURCE_DAMAGE_ISODUCTILE_ID) then
source_damage_isoDuctile_offset(p) = sourceOffset
exit
endif
enddo
if (all(phase_source(:,p) /= SOURCE_DAMAGE_ISODUCTILE_ID)) cycle
associate(prm => param(source_damage_isoDuctile_instance(p)), &
config => config_phase(p))
prm%aTol = config%getFloat('isoductile_atol',defaultVal = 1.0e-3_pReal)
prm%N = config%getFloat('isoductile_ratesensitivity')
prm%critPlasticStrain = config%getFloat('isoductile_criticalplasticstrain')
! sanity checks
if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isoductile_atol'
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_ratesensitivity'
if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain'
@ -104,33 +89,18 @@ subroutine source_damage_isoDuctile_init
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
case ('isoductile_drivingforce')
prm%outputID = [prm%outputID, damage_drivingforce_ID]
end select
enddo
NofMyPhase=count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
end associate
phase = p
NofMyPhase=count(material_phaseAt==phase) * discretization_nIP
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
enddo
end subroutine source_damage_isoDuctile_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
@ -140,23 +110,28 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: &
phase, constituent, instance, homog, sourceOffset, damageOffset
phase, &
constituent, &
sourceOffset, &
damageOffset, &
homog
phase = material_phaseAt(ipc,el)
constituent = material_phasememberAt(ipc,ip,el)
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
associate(prm => param(source_damage_isoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sum(plasticState(phase)%slipRate(:,constituent))/ &
((damage(homog)%p(damageOffset))**param(instance)%N)/ &
param(instance)%critPlasticStrain
sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain
end associate
end subroutine source_damage_isoDuctile_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
@ -170,16 +145,17 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
sourceOffset
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)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
end subroutine source_damage_isoDuctile_getRateAndItsTangent
@ -190,15 +166,14 @@ subroutine source_damage_isoDuctile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
integer :: sourceOffset, o, instance
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
integer :: o
associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_drivingforce_ID)
associate(prm => param(source_damage_isoDuctile_instance(phase)), &
stt => sourceState(phase)%p(source_damage_isoDuctile_offset(phase))%state)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('isoductile_drivingforce')
call results_writeDataset(group,stt,'tbd','driving force','tbd')
end select
enddo outputsLoop
@ -206,5 +181,4 @@ subroutine source_damage_isoDuctile_results(phase,group)
end subroutine source_damage_isoDuctile_results
end module source_damage_isoDuctile

55
src/source_thermal_dissipation.f90
View File

@ -39,46 +39,46 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_dissipation_init
integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p
integer :: Ninstance,sourceOffset,NofMyPhase,p
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'; flush(6)
Ninstance = count(phase_source == SOURCE_thermal_dissipation_ID)
if (Ninstance == 0) return
Ninstance = count(phase_source == SOURCE_THERMAL_DISSIPATION_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_thermal_dissipation_offset(material_Nphase), source=0)
allocate(source_thermal_dissipation_instance(material_Nphase), source=0)
allocate(source_thermal_dissipation_offset (size(config_phase)), source=0)
allocate(source_thermal_dissipation_instance(size(config_phase)), source=0)
allocate(param(Ninstance))
do p = 1, material_Nphase
source_thermal_dissipation_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_dissipation_ID)
do source = 1, phase_Nsources(p)
if (phase_source(source,p) == SOURCE_thermal_dissipation_ID) &
source_thermal_dissipation_offset(p) = source
enddo
enddo
do p = 1, size(config_phase)
if (all(phase_source(:,p) /= SOURCE_THERMAL_DISSIPATION_ID)) cycle
instance = source_thermal_dissipation_instance(p)
param(instance)%kappa = config_phase(p)%getFloat('dissipation_coldworkcoeff')
NofMyPhase = count(material_phaseAt==p) * discretization_nIP
sourceOffset = source_thermal_dissipation_offset(p)
source_thermal_dissipation_instance(p) = count(phase_source(:,1:p) == SOURCE_THERMAL_DISSIPATION_ID)
do sourceOffset = 1, phase_Nsources(p)
if (phase_source(sourceOffset,p) == SOURCE_THERMAL_DISSIPATION_ID) then
source_thermal_dissipation_offset(p) = sourceOffset
exit
endif
enddo
if (all(phase_source(:,p) /= SOURCE_THERMAL_DISSIPATION_ID)) cycle
associate(prm => param(source_thermal_dissipation_instance(p)), &
config => config_phase(p))
prm%kappa = config%getFloat('dissipation_coldworkcoeff')
NofMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NofMyPhase,0,0,0)
end associate
enddo
end subroutine source_thermal_dissipation_init
!--------------------------------------------------------------------------------------------------
!> @brief returns dissipation rate
!> @brief Ninstances dissipation rate
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase)
subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase)
integer, intent(in) :: &
phase
@ -86,16 +86,15 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar
Tstar
real(pReal), intent(in), dimension(3,3) :: &
Lp
real(pReal), intent(out) :: &
TDot, &
dTDOT_dT
integer :: &
instance
dTDot_dT
instance = source_thermal_dissipation_instance(phase)
TDot = param(instance)%kappa*sum(abs(Tstar*Lp))
dTDOT_dT = 0.0_pReal
associate(prm => param(source_thermal_dissipation_instance(phase)))
TDot = prm%kappa*sum(abs(Tstar*Lp))
dTDot_dT = 0.0_pReal
end associate
end subroutine source_thermal_dissipation_getRateAndItsTangent

68
src/source_thermal_externalheat.f90
View File

@ -43,43 +43,40 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_init
integer :: maxNinstance,instance,source,sourceOffset,NofMyPhase,p
integer :: Ninstance,sourceOffset,NofMyPhase,p
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'; flush(6)
maxNinstance = count(phase_source == SOURCE_thermal_externalheat_ID)
if (maxNinstance == 0) return
Ninstance = count(phase_source == SOURCE_thermal_externalheat_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_thermal_externalheat_offset(material_Nphase), source=0)
allocate(source_thermal_externalheat_instance(material_Nphase), source=0)
do p = 1, material_Nphase
source_thermal_externalheat_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_externalheat_ID)
do source = 1, phase_Nsources(p)
if (phase_source(source,p) == SOURCE_thermal_externalheat_ID) &
source_thermal_externalheat_offset(p) = source
enddo
enddo
allocate(param(maxNinstance))
allocate(source_thermal_externalheat_offset (size(config_phase)), source=0)
allocate(source_thermal_externalheat_instance(size(config_phase)), source=0)
allocate(param(Ninstance))
do p = 1, size(config_phase)
source_thermal_externalheat_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_externalheat_ID)
do sourceOffset = 1, phase_Nsources(p)
if (phase_source(sourceOffset,p) == SOURCE_thermal_externalheat_ID) then
source_thermal_externalheat_offset(p) = sourceOffset
exit
endif
enddo
if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle
instance = source_thermal_externalheat_instance(p)
sourceOffset = source_thermal_externalheat_offset(p)
associate(prm => param(source_thermal_externalheat_instance(p)), &
config => config_phase(p))
prm%time = config%getFloats('externalheat_time')
prm%nIntervals = size(prm%time) - 1
prm%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(prm%time))
NofMyPhase = count(material_phaseAt==p) * discretization_nIP
param(instance)%time = config_phase(p)%getFloats('externalheat_time')
param(instance)%nIntervals = size(param(instance)%time) - 1
param(instance)%heat_rate = config_phase(p)%getFloats('externalheat_rate',requiredSize = size(param(instance)%time))
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
end associate
enddo
end subroutine source_thermal_externalheat_init
@ -94,6 +91,7 @@ subroutine source_thermal_externalheat_dotState(phase, of)
integer, intent(in) :: &
phase, &
of
integer :: &
sourceOffset
@ -115,25 +113,27 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas
real(pReal), intent(out) :: &
TDot, &
dTDot_dT
integer :: &
instance, sourceOffset, interval
sourceOffset, interval
real(pReal) :: &
frac_time
instance = source_thermal_externalheat_instance(phase)
sourceOffset = source_thermal_externalheat_offset(phase)
do interval = 1, param(instance)%nIntervals ! scan through all rate segments
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - param(instance)%time(interval)) &
/ (param(instance)%time(interval+1) - param(instance)%time(interval)) ! fractional time within segment
associate(prm => param(source_thermal_externalheat_instance(phase)))
do interval = 1, prm%nIntervals ! scan through all rate segments
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%time(interval)) &
/ (prm%time(interval+1) - prm%time(interval)) ! fractional time within segment
if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == param(instance)%nIntervals) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = param(instance)%heat_rate(interval ) * (1.0_pReal - frac_time) + &
param(instance)%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
TDot = prm%heat_rate(interval ) * (1.0_pReal - frac_time) + &
prm%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside of bounds
enddo
dTDot_dT = 0.0
end associate
end subroutine source_thermal_externalheat_getRateAndItsTangent

27
src/thermal_conduction.f90
View File

@ -26,7 +26,7 @@ module thermal_conduction
public :: &
thermal_conduction_init, &
thermal_conduction_getSourceAndItsTangent, &
thermal_conduction_getConductivity33, &
thermal_conduction_getConductivity, &
thermal_conduction_getSpecificHeat, &
thermal_conduction_getMassDensity, &
thermal_conduction_putTemperatureAndItsRate, &
@ -41,17 +41,14 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_init
integer :: maxNinstance,NofMyHomog,h
integer :: Ninstance,NofMyHomog,h
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'; flush(6)
maxNinstance = count(thermal_type == THERMAL_conduction_ID)
if (maxNinstance == 0) return
Ninstance = count(thermal_type == THERMAL_conduction_ID)
allocate(param(Ninstance))
allocate(param(maxNinstance))
do h = 1, size(thermal_type)
do h = 1, size(config_homogenization)
if (thermal_type(h) /= THERMAL_conduction_ID) cycle
associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h))
@ -137,27 +134,27 @@ end subroutine thermal_conduction_getSourceAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized thermal conductivity in reference configuration
!--------------------------------------------------------------------------------------------------
function thermal_conduction_getConductivity33(ip,el)
function thermal_conduction_getConductivity(ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
thermal_conduction_getConductivity33
thermal_conduction_getConductivity
integer :: &
grain
thermal_conduction_getConductivity33 = 0.0_pReal
thermal_conduction_getConductivity = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_conduction_getConductivity33 = thermal_conduction_getConductivity33 + &
crystallite_push33ToRef(grain,ip,el,lattice_thermalConductivity33(:,:,material_phaseAt(grain,el)))
thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
crystallite_push33ToRef(grain,ip,el,lattice_thermalConductivity(:,:,material_phaseAt(grain,el)))
enddo
thermal_conduction_getConductivity33 = thermal_conduction_getConductivity33 &
thermal_conduction_getConductivity = thermal_conduction_getConductivity &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getConductivity33
end function thermal_conduction_getConductivity
!--------------------------------------------------------------------------------------------------

5
src/thermal_isothermal.f90
View File

@ -7,10 +7,7 @@ module thermal_isothermal
use material
implicit none
private
public :: &
thermal_isothermal_init
public
contains