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

Merge branch 'development' into HDF5-results

This commit is contained in:
Martin Diehl 2020-03-13 23:32:25 +01:00
parent 81e98055dd c9754767e4
commit d707378820
29 changed files with 1350 additions and 1692 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
VERSION
View File

@ -1 +1 @@
v2.0.3-1862-g0b340a6d v2.0.3-1920-ge8c16e7f

2
examples/ConfigFiles/Phase_Isotropic_AluminumIsotropic.config
View File

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

2
examples/ConfigFiles/Phase_Isotropic_FreeSurface.config
View File

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

3
examples/ConfigFiles/Phase_None_IsotropicVolumePreservation.config
View File

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

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

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

21
src/CPFEM.f90
View File

@ -27,9 +27,6 @@ module CPFEM
implicit none implicit none
private 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 :: & real(pReal), dimension (:,:,:), allocatable, private :: &
CPFEM_cs !< Cauchy stress CPFEM_cs !< Cauchy stress
real(pReal), dimension (:,:,:,:), allocatable, private :: & 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 integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog, mySource 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) 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 elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress CPFEM_cs(1:6,ip,elCP) = rnd * ODD_STRESS
CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
chosenThermal2: select case (thermal_type(material_homogenizationAt(elCP))) chosenThermal2: select case (thermal_type(material_homogenizationAt(elCP)))
case (THERMAL_conduction_ID) chosenThermal2 case (THERMAL_conduction_ID) chosenThermal2
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & 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 endif
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = rnd*CPFEM_odd_stress CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian*math_identity2nd(6) CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
!*** deformation gradient is not outdated !*** deformation gradient is not outdated
@ -257,8 +257,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
else terminalIllness else terminalIllness
@ -331,6 +331,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
end subroutine CPFEM_general end subroutine CPFEM_general
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Forward data for new time increment. !> @brief Forward data for new time increment.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

2
src/DAMASK_marc.f90
View File

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

14
src/HDF5_utilities.f90
View File

@ -77,20 +77,6 @@ module HDF5_utilities
module procedure HDF5_addAttribute_real_array module procedure HDF5_addAttribute_real_array
end interface HDF5_addAttribute 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 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'),& prm%n0_sl = lattice_slip_normal(prm%N_sl,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) 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) .and. (prm%N_sl(1) == 12)
if(prm%fccTwinTransNucleation) & 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_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl))
prm%rho_dip_0 = config%getFloats('rhoedgedip0',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) & prm%omega = config%getFloat('omega', defaultVal = 1000.0_pReal) &
* merge(12.0_pReal, & * merge(12.0_pReal, &
8.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 ! expand: family => system
@ -335,7 +335,7 @@ module subroutine plastic_dislotwin_init
config%getFloat('a_bcc', defaultVal=0.0_pReal), & config%getFloat('a_bcc', defaultVal=0.0_pReal), &
config%getFloat('a_fcc', 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 = config%getFloats('ndot0_trans')
prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr) prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr)
endif endif

42
src/damage_local.f90
View File

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

60
src/damage_nonlocal.f90
View File

@ -18,15 +18,9 @@ module damage_nonlocal
implicit none implicit none
private private
enum, bind(c)
enumerator :: &
undefined_ID, &
damage_ID
end enum
type :: tParameters type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: & character(len=pStringLen), allocatable, dimension(:) :: &
outputID output
end type tParameters end type tParameters
type(tparameters), dimension(:), allocatable :: & type(tparameters), dimension(:), allocatable :: &
@ -35,10 +29,10 @@ module damage_nonlocal
public :: & public :: &
damage_nonlocal_init, & damage_nonlocal_init, &
damage_nonlocal_getSourceAndItsTangent, & damage_nonlocal_getSourceAndItsTangent, &
damage_nonlocal_getDiffusion33, & damage_nonlocal_getDiffusion, &
damage_nonlocal_getMobility, & damage_nonlocal_getMobility, &
damage_nonlocal_putNonLocalDamage, & damage_nonlocal_putNonLocalDamage, &
damage_nonlocal_Results damage_nonlocal_results
contains contains
@ -48,29 +42,18 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_init subroutine damage_nonlocal_init
integer :: maxNinstance,o,NofMyHomog,h integer :: Ninstance,NofMyHomog,h
character(len=pStringLen), dimension(:), allocatable :: outputs
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'; flush(6)
maxNinstance = count(damage_type == DAMAGE_nonlocal_ID) Ninstance = count(damage_type == DAMAGE_nonlocal_ID)
if (maxNinstance == 0) return allocate(param(Ninstance))
allocate(param(maxNinstance)) do h = 1, size(config_homogenization)
do h = 1, size(damage_type)
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
associate(prm => param(damage_typeInstance(h)),config => config_homogenization(h)) associate(prm => param(damage_typeInstance(h)),config => config_homogenization(h))
outputs = config%getStrings('(output)',defaultVal=emptyStringArray) prm%output = 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
NofMyHomog = count(material_homogenizationAt == h) NofMyHomog = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1 damageState(h)%sizeState = 1
@ -145,28 +128,28 @@ end subroutine damage_nonlocal_getSourceAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized non local damage diffusion tensor in reference configuration !> @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) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
damage_nonlocal_getDiffusion33 damage_nonlocal_getDiffusion
integer :: & integer :: &
homog, & homog, &
grain grain
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damage_nonlocal_getDiffusion33 = 0.0_pReal damage_nonlocal_getDiffusion = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog) do grain = 1, homogenization_Ngrains(homog)
damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + & damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phaseAt(grain,el))) crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion(1:3,1:3,material_phaseAt(grain,el)))
enddo enddo
damage_nonlocal_getDiffusion33 = & damage_nonlocal_getDiffusion = &
charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal) 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 integer, intent(in) :: homog
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer :: o integer :: o
associate(prm => param(damage_typeInstance(homog))) associate(prm => param(damage_typeInstance(homog)))
outputsLoop: do o = 1,size(prm%output)
outputsLoop: do o = 1,size(prm%outputID) select case(prm%output(o))
select case(prm%outputID(o)) case ('damage')
case (damage_ID)
call results_writeDataset(group,damage(homog)%p,'phi',& call results_writeDataset(group,damage(homog)%p,'phi',&
'damage indicator','-') 'damage indicator','-')
end select 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 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 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)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward call utilities_FFTvectorForward
@ -294,7 +294,7 @@ subroutine updateReference
mu_ref = 0.0_pReal mu_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 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) mu_ref = mu_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo enddo; enddo; enddo
K_ref = K_ref*wgt 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 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 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)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward call utilities_FFTvectorForward
@ -294,9 +294,8 @@ subroutine updateReference
mu_ref = 0.0_pReal mu_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
K_ref = K_ref + thermal_conduction_getConductivity33(1,cell) K_ref = K_ref + thermal_conduction_getConductivity(1,cell)
mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* & mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* thermal_conduction_getSpecificHeat(1,cell)
thermal_conduction_getSpecificHeat(1,cell)
enddo; enddo; enddo enddo; enddo; enddo
K_ref = K_ref*wgt K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) 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, & sdot0, &
n n
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
critDisp, &
critLoad critLoad
end type real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems
end type tParameters
! Begin Deprecated type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
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
public :: & public :: &
kinematics_cleavage_opening_init, & kinematics_cleavage_opening_init, &
@ -60,61 +47,54 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_init subroutine kinematics_cleavage_opening_init
integer, allocatable, dimension(:) :: tempInt integer :: Ninstance,p
real(pReal), allocatable, dimension(:) :: tempFloat character(len=pStringLen) :: extmsg = ''
integer :: maxNinstance,p,instance
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'; flush(6)
maxNinstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID) Ninstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID)
if (maxNinstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0) allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0)
do p = 1, size(config_phase) allocate(param(Ninstance))
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)
do p = 1, size(config_phase) 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 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)) associate(prm => param(kinematics_cleavage_opening_instance(p)), &
kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat config => config_phase(p))
tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt)) prm%Ncleavage = config%getInts('ncleavage')
kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat prm%totalNcleavage = sum(prm%Ncleavage)
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = & prm%n = config%getFloat('anisobrittle_ratesensitivity')
min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,p),& ! limit active cleavage systems per family to min of available and requested prm%sdot0 = config%getFloat('anisobrittle_sdot0')
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance))
kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether prm%critLoad = config%getFloats('anisobrittle_criticalload',requiredSize=size(prm%Ncleavage))
if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')') prm%cleavage_systems = lattice_SchmidMatrix_cleavage(prm%Ncleavage,config%getString('lattice_structure'),&
if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) & config%getFloat('c/a',defaultVal=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)) & ! expand: family => system
call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')') prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage)
if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')') ! 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 enddo
end subroutine kinematics_cleavage_opening_init end subroutine kinematics_cleavage_opening_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -130,72 +110,54 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld !< damage velocity gradient Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
integer :: & integer :: &
instance, phase, &
homog, damageOffset, & homog, damageOffset, &
f, i, index_myFamily, k, l, m, n i, k, l, m, n
real(pReal) :: & real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit, & traction_d, traction_t, traction_n, traction_crit, &
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
phase = material_phaseAt(ipc,el)
instance = kinematics_cleavage_opening_instance(phase)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el) damageOffset = damageMapping(homog)%p(ip,el)
Ld = 0.0_pReal Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal dLd_dTstar = 0.0_pReal
do f = 1,lattice_maxNcleavageFamily associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el))))
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family do i = 1,prm%totalNcleavage
do i = 1,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset)**2.0_pReal
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_d = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,1,i))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) if (abs(traction_d) > traction_crit + tol_math_check) then
traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* & udotd = sign(1.0_pReal,traction_d)* prm%sdot0 * ((abs(traction_d) - traction_crit)/traction_crit)**prm%n
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i)
udotd = & dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%n / (abs(traction_d) - traction_crit)
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:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* & + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i)
lattice_Scleavage(m,n,1,index_myFamily+i,phase)
endif endif
udott = & traction_t = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,2,i))
sign(1.0_pReal,traction_t)* & if (abs(traction_t) > traction_crit + tol_math_check) then
kinematics_cleavage_opening_sdot_0(instance)* & udott = sign(1.0_pReal,traction_t)* prm%sdot0 * ((abs(traction_t) - traction_crit)/traction_crit)**prm%n
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i)
if (abs(udott) > tol_math_check) then dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%n / (abs(traction_t) - traction_crit)
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:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* & + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i)
lattice_Scleavage(m,n,2,index_myFamily+i,phase)
endif endif
udotn = & traction_n = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,3,i))
sign(1.0_pReal,traction_n)* & if (abs(traction_n) > traction_crit + tol_math_check) then
kinematics_cleavage_opening_sdot_0(instance)* & udotn = sign(1.0_pReal,traction_n)* prm%sdot0 * ((abs(traction_n) - traction_crit)/traction_crit)**prm%n
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i)
if (abs(udotn) > tol_math_check) then dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%n / (abs(traction_n) - traction_crit)
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:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* & + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i)
lattice_Scleavage(m,n,3,index_myFamily+i,phase)
endif endif
enddo enddo
enddo end associate
end subroutine kinematics_cleavage_opening_LiAndItsTangent end subroutine kinematics_cleavage_opening_LiAndItsTangent

107
src/kinematics_slipplane_opening.f90
View File

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

85
src/kinematics_thermal_expansion.f90
View File

@ -15,9 +15,13 @@ module kinematics_thermal_expansion
implicit none implicit none
private private
integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance
type :: tParameters type :: tParameters
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal) :: &
expansion T_ref
real(pReal), dimension(3,3,3) :: &
expansion = 0.0_pReal
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param type(tParameters), dimension(:), allocatable :: param
@ -36,35 +40,39 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_init subroutine kinematics_thermal_expansion_init
integer :: & integer :: Ninstance,p,i
Ninstance, & real(pReal), dimension(:), allocatable :: temp
p, i
real(pReal), dimension(:), allocatable :: &
temp
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'; flush(6)
Ninstance = count(phase_kinematics == KINEMATICS_thermal_expansion_ID) Ninstance = count(phase_kinematics == KINEMATICS_thermal_expansion_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(kinematics_thermal_expansion_instance(size(config_phase)), source=0)
allocate(param(Ninstance)) 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 if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle
! ToDo: Here we need to decide how to extend the concept of instances to associate(prm => param(kinematics_thermal_expansion_instance(p)), &
! kinetics and sources. I would suggest that the same mechanism exists at maximum once per phase 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) ! read up to three parameters (constant, linear, quadratic with T)
temp = config_phase(p)%getFloats('thermal_expansion11') temp = config%getFloats('thermal_expansion11')
!lattice_thermalExpansion33(1,1,1:size(temp),p) = temp prm%expansion(1,1,1:size(temp)) = temp
temp = config_phase(p)%getFloats('thermal_expansion22', & temp = config%getFloats('thermal_expansion22',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) prm%expansion(2,2,1:size(temp)) = temp
!lattice_thermalExpansion33(2,2,1:size(temp),p) = temp temp = config%getFloats('thermal_expansion33',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
temp = config_phase(p)%getFloats('thermal_expansion33', & prm%expansion(3,3,1:size(temp)) = temp
defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(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 enddo
end subroutine kinematics_thermal_expansion_init end subroutine kinematics_thermal_expansion_init
@ -77,18 +85,18 @@ pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset)
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &
homog, offset homog, &
offset
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though)
associate(prm => param(kinematics_thermal_expansion_instance(phase)))
kinematics_thermal_expansion_initialStrain = & kinematics_thermal_expansion_initialStrain = &
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * & (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient
lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * & (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient
lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient end associate
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * &
lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient
end function kinematics_thermal_expansion_initialStrain end function kinematics_thermal_expansion_initialStrain
@ -106,29 +114,30 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip,
Li !< thermal velocity gradient Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & 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) dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
integer :: & integer :: &
phase, & phase, &
homog, offset homog
real(pReal) :: & real(pReal) :: &
T, TRef, TDot T, TDot
phase = material_phaseAt(ipc,el) phase = material_phaseAt(ipc,el)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
offset = thermalMapping(homog)%p(ip,el) T = temperature(homog)%p(thermalMapping(homog)%p(ip,el))
T = temperature(homog)%p(offset) TDot = temperatureRate(homog)%p(thermalMapping(homog)%p(ip,el))
TDot = temperatureRate(homog)%p(offset)
TRef = lattice_referenceTemperature(phase)
associate(prm => param(kinematics_thermal_expansion_instance(phase)))
Li = TDot * ( & Li = TDot * ( &
lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient + prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient + prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
) / & ) / &
(1.0_pReal & (1.0_pReal &
+ lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. & + prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & + prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
) )
end associate
dLi_dTstar = 0.0_pReal dLi_dTstar = 0.0_pReal
end subroutine kinematics_thermal_expansion_LiAndItsTangent 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 implicit none
public public
! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds ! 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) #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 #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 #endif
integer, parameter, public :: pLongInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) integer, parameter :: 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 :: pStringLen = 256 !< default string length
integer, parameter, public :: pPathLen = 4096 !< maximum length of a path name on linux 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 real(pReal), dimension(:), pointer :: p
end type group_float end type group_float
type, public :: group_int type :: group_int
integer, dimension(:), pointer :: p integer, dimension(:), pointer :: p
end type group_int end type group_int
! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array ! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array
type, public :: tState type :: tState
integer :: & integer :: &
sizeState = 0, & !< size of state sizeState = 0, & !< size of state
sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
@ -57,7 +58,7 @@ module prec
RKCK45dotState RKCK45dotState
end type end type
type, extends(tState), public :: tPlasticState type, extends(tState) :: tPlasticState
logical :: & logical :: &
nonlocal = .false. nonlocal = .false.
real(pReal), pointer, dimension(:,:) :: & real(pReal), pointer, dimension(:,:) :: &
@ -65,22 +66,22 @@ module prec
accumulatedSlip !< accumulated plastic slip accumulatedSlip !< accumulated plastic slip
end type end type
type, public :: tSourceState type :: tSourceState
type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase
end type end type
type, public :: tHomogMapping type :: tHomogMapping
integer, pointer, dimension(:,:) :: p integer, pointer, dimension(:,:) :: p
end type 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_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 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::] emptyIntArray = [integer::]
real(pReal), dimension(0), parameter, public :: & real(pReal), dimension(0), parameter :: &
emptyRealArray = [real(pReal)::] emptyRealArray = [real(pReal)::]
character(len=pStringLen), dimension(0), parameter, public :: & character(len=pStringLen), dimension(0), parameter :: &
emptyStringArray = [character(len=pStringLen)::] emptyStringArray = [character(len=pStringLen)::]
private :: & private :: &

11
src/quaternions.f90
View File

@ -10,7 +10,7 @@ module quaternions
use IO use IO
implicit none implicit none
public private
real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion. real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion.
@ -96,8 +96,12 @@ module quaternions
module procedure aimag__ module procedure aimag__
end interface aimag end interface aimag
private :: & public :: &
unitTest quaternions_init, &
assignment(=), &
conjg, aimag, &
log, exp, &
real
contains contains
@ -511,6 +515,7 @@ subroutine unitTest
q_2 = conjg(q_2) - inverse(q_2) 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') if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(0,ext_msg='inverse/conjg')
endif 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 !(defined(__GFORTRAN__) && __GNUC__ < 9)
if (norm2(aimag(q)) > 0.0_pReal) then 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 debug
use IO use IO
use math use math
use material
use discretization use discretization
use material
use config use config
use lattice use lattice
use results use results
@ -22,15 +22,6 @@ module source_damage_anisoBrittle
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
source_damage_anisoBrittle_instance !< instance of source mechanism source_damage_anisoBrittle_instance !< instance of source mechanism
integer, 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 type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
aTol, & aTol, &
@ -45,8 +36,8 @@ module source_damage_anisoBrittle
totalNcleavage totalNcleavage
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
Ncleavage Ncleavage
integer(kind(undefined_ID)), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
outputID !< ID of each post result output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
@ -67,56 +58,39 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_init subroutine source_damage_anisoBrittle_init
integer :: Ninstance,phase,instance,source,sourceOffset integer :: Ninstance,sourceOffset,NofMyPhase,p
integer :: NofMyPhase,p ,i character(len=pStringLen) :: extmsg = ''
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'; flush(6)
Ninstance = count(phase_source == SOURCE_damage_anisoBrittle_ID) Ninstance = count(phase_source == SOURCE_DAMAGE_ANISOBRITTLE_ID)
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0) allocate(source_damage_anisoBrittle_offset (size(config_phase)), source=0)
allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0) allocate(source_damage_anisoBrittle_instance(size(config_phase)), 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(param(Ninstance)) allocate(param(Ninstance))
do p = 1, size(config_phase) 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 if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISOBRITTLE_ID)) cycle
associate(prm => param(source_damage_anisoBrittle_instance(p)), & associate(prm => param(source_damage_anisoBrittle_instance(p)), &
config => config_phase(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%N = config%getFloat('anisobrittle_ratesensitivity')
prm%sdot_0 = config%getFloat('anisobrittle_sdot0') 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%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage))
prm%critLoad = config%getFloats('anisobrittle_criticalload', 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%critDisp = math_expand(prm%critDisp, prm%Ncleavage)
prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage)
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload' ! sanity checks
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement' 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 ! exit if any parameter is out of range
if (extmsg /= '') & if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output pararameters ! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray) prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('anisobrittle_drivingforce') NofMyPhase = count(material_phaseAt==p) * discretization_nIP
prm%outputID = [prm%outputID, damage_drivingforce_ID] call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
end select
enddo
end associate 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 enddo
end subroutine source_damage_anisoBrittle_init end subroutine source_damage_anisoBrittle_init
@ -178,49 +137,42 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
el !< element el !< element
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S S
integer :: & integer :: &
phase, & phase, &
constituent, & constituent, &
instance, &
sourceOffset, & sourceOffset, &
damageOffset, & damageOffset, &
homog, & homog, &
f, i, index_myFamily, index i
real(pReal) :: & real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit traction_d, traction_t, traction_n, traction_crit
phase = material_phaseAt(ipc,el) phase = material_phaseAt(ipc,el)
constituent = material_phasememberAt(ipc,ip,el) constituent = material_phasememberAt(ipc,ip,el)
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase) sourceOffset = source_damage_anisoBrittle_offset(phase)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,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 sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
do i = 1, prm%totalNcleavage
index = 1 traction_d = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,1,i))
do f = 1,lattice_maxNcleavageFamily traction_t = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,2,i))
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family traction_n = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,3,i))
do i = 1,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) traction_crit = prm%critLoad(i)*damage(homog)%p(damageOffset)**2.0_pReal
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 = param(instance)%critLoad(index)* & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) = 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 enddo
end associate
end subroutine source_damage_anisoBrittle_dotState end subroutine source_damage_anisoBrittle_dotState
@ -238,16 +190,17 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer :: & integer :: &
sourceOffset sourceOffset
sourceOffset = source_damage_anisoBrittle_offset(phase) 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) dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
end subroutine source_damage_anisoBrittle_getRateAndItsTangent end subroutine source_damage_anisoBrittle_getRateAndItsTangent
@ -258,15 +211,14 @@ subroutine source_damage_anisoBrittle_results(phase,group)
integer, intent(in) :: phase integer, intent(in) :: phase
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer :: sourceOffset, o, instance
instance = source_damage_anisoBrittle_instance(phase) integer :: o
sourceOffset = source_damage_anisoBrittle_offset(phase)
associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state) associate(prm => param(source_damage_anisoBrittle_instance(phase)), &
outputsLoop: do o = 1,size(prm%outputID) stt => sourceState(phase)%p(source_damage_anisoBrittle_offset(phase))%state)
select case(prm%outputID(o)) outputsLoop: do o = 1,size(prm%output)
case (damage_drivingforce_ID) select case(trim(prm%output(o)))
case ('anisobrittle_drivingforce')
call results_writeDataset(group,stt,'tbd','driving force','tbd') call results_writeDataset(group,stt,'tbd','driving force','tbd')
end select end select
enddo outputsLoop 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_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism source_damage_anisoDuctile_instance !< instance of damage source mechanism
enum, bind(c)
enumerator :: undefined_ID, &
damage_drivingforce_ID
end enum
type, private :: tParameters !< container type for internal constitutive parameters type, private :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
aTol, & aTol, &
@ -37,8 +31,8 @@ module source_damage_anisoDuctile
totalNslip totalNslip
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
Nslip Nslip
integer(kind(undefined_ID)), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
outputID output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
@ -59,58 +53,46 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_init subroutine source_damage_anisoDuctile_init
integer :: Ninstance,phase,instance,source,sourceOffset integer :: Ninstance,sourceOffset,NofMyPhase,p
integer :: NofMyPhase,p ,i character(len=pStringLen) :: extmsg = ''
integer(kind(undefined_ID)) :: &
outputID
character(len=pStringLen) :: &
extmsg = ''
character(len=pStringLen), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'; flush(6)
Ninstance = count(phase_source == SOURCE_damage_anisoDuctile_ID) Ninstance = count(phase_source == SOURCE_DAMAGE_ANISODUCTILE_ID)
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoDuctile_offset (size(config_phase)), source=0) allocate(source_damage_anisoDuctile_offset (size(config_phase)), source=0)
allocate(source_damage_anisoDuctile_instance(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)) allocate(param(Ninstance))
do p = 1, size(config_phase) 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 if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISODUCTILE_ID)) cycle
associate(prm => param(source_damage_anisoDuctile_instance(p)), & associate(prm => param(source_damage_anisoDuctile_instance(p)), &
config => config_phase(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%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)) prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip))
! expand: family => system ! expand: family => system
prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip) 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' if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -119,30 +101,13 @@ subroutine source_damage_anisoDuctile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output pararameters ! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray) prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('anisoductile_drivingforce') NofMyPhase=count(material_phaseAt==p) * discretization_nIP
prm%outputID = [prm%outputID, damage_drivingforce_ID] call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
end select
enddo
end associate 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 enddo
end subroutine source_damage_anisoDuctile_init end subroutine source_damage_anisoDuctile_init
@ -157,28 +122,28 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
integer :: & integer :: &
phase, & phase, &
constituent, & constituent, &
sourceOffset, & sourceOffset, &
homog, damageOffset, & damageOffset, &
instance, & homog, &
i i
phase = material_phaseAt(ipc,el) phase = material_phaseAt(ipc,el)
constituent = material_phasememberAt(ipc,ip,el) constituent = material_phasememberAt(ipc,ip,el)
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase) sourceOffset = source_damage_anisoDuctile_offset(phase)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el) damageOffset = damageMapping(homog)%p(ip,el)
associate(prm => param(source_damage_anisoDuctile_instance(phase)))
do i = 1, param(instance)%totalNslip do i = 1, prm%totalNslip
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + & = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
plasticState(phase)%slipRate(i,constituent)/ & + plasticState(phase)%slipRate(i,constituent)/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain(i)
((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i)
enddo enddo
end associate
end subroutine source_damage_anisoDuctile_dotState end subroutine source_damage_anisoDuctile_dotState
@ -196,16 +161,17 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer :: & integer :: &
sourceOffset sourceOffset
sourceOffset = source_damage_anisoDuctile_offset(phase) 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) dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
end subroutine source_damage_anisoDuctile_getRateAndItsTangent end subroutine source_damage_anisoDuctile_getRateAndItsTangent
@ -216,15 +182,14 @@ subroutine source_damage_anisoDuctile_results(phase,group)
integer, intent(in) :: phase integer, intent(in) :: phase
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer :: sourceOffset, o, instance
instance = source_damage_anisoDuctile_instance(phase) integer :: o
sourceOffset = source_damage_anisoDuctile_offset(phase)
associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state) associate(prm => param(source_damage_anisoDuctile_instance(phase)), &
outputsLoop: do o = 1,size(prm%outputID) stt => sourceState(phase)%p(source_damage_anisoDuctile_offset(phase))%state)
select case(prm%outputID(o)) outputsLoop: do o = 1,size(prm%output)
case (damage_drivingforce_ID) select case(trim(prm%output(o)))
case ('anisoductile_drivingforce')
call results_writeDataset(group,stt,'tbd','driving force','tbd') call results_writeDataset(group,stt,'tbd','driving force','tbd')
end select end select
enddo outputsLoop enddo outputsLoop

121
src/source_damage_isoBrittle.f90
View File

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

110
src/source_damage_isoDuctile.f90
View File

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

55
src/source_thermal_dissipation.f90
View File

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

68
src/source_thermal_externalheat.f90
View File

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

27
src/thermal_conduction.f90
View File

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