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

polishing: adding pInt, removing unused use-statements etc

This commit is contained in:
Martin Diehl 2012-02-21 16:00:00 +00:00
parent 23cda48709
commit 79663a7f76
5 changed files with 544 additions and 573 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

85
code/constitutive.f90
View File

@ -68,8 +68,9 @@ CONTAINS
!* Module initialization *
!**************************************
subroutine constitutive_init()
use prec, only: pReal,pInt
use debug, only: debug_verbosity, debug_selectiveDebugger, debug_e, debug_i, debug_g
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt
use debug, only: debug_verbosity
use numerics, only: numerics_integrator
use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
@ -82,7 +83,7 @@ use constitutive_nonlocal
implicit none
integer(pInt), parameter :: fileunit = 200
integer(pInt), parameter :: fileunit = 200_pInt
integer(pInt) g, & ! grain number
i, & ! integration point number
e, & ! element number
@ -114,7 +115,7 @@ close(fileunit)
! --- WRITE DESCRIPTION FILE FOR CONSTITUTIVE PHASE OUTPUT ---
call IO_write_jobFile(fileunit,'outputConstitutive')
do p = 1,material_Nphase
do p = 1_pInt,material_Nphase
i = phase_constitutionInstance(p) ! which instance of a constitution is present phase
knownConstitution = .true. ! assume valid
select case(phase_constitution(p)) ! split per constitiution
@ -141,7 +142,7 @@ do p = 1,material_Nphase
write(fileunit,*)
if (knownConstitution) then
write(fileunit,'(a)') '(constitution)'//char(9)//trim(phase_constitution(p))
do e = 1,phase_Noutput(p)
do e = 1_pInt,phase_Noutput(p)
write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
@ -166,22 +167,22 @@ allocate(constitutive_aTolState(gMax,iMax,eMax))
allocate(constitutive_sizeDotState(gMax,iMax,eMax)) ; constitutive_sizeDotState = 0_pInt
allocate(constitutive_sizeState(gMax,iMax,eMax)) ; constitutive_sizeState = 0_pInt
allocate(constitutive_sizePostResults(gMax,iMax,eMax)); constitutive_sizePostResults = 0_pInt
if (any(numerics_integrator == 1)) then
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(gMax,iMax,eMax))
allocate(constitutive_previousDotState2(gMax,iMax,eMax))
endif
if (any(numerics_integrator == 4)) then
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(gMax,iMax,eMax))
endif
if (any(numerics_integrator == 5)) then
if (any(numerics_integrator == 5_pInt)) then
allocate(constitutive_RKCK45dotState(6,gMax,iMax,eMax))
endif
!$OMP PARALLEL DO PRIVATE(myNgrains,myInstance)
do e = 1,mesh_NcpElems ! loop over elements
do e = 1_pInt,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
do g = 1,myNgrains ! loop over grains
do i = 1_pInt,FE_Nips(mesh_element(2,e)) ! loop over IPs
do g = 1_pInt,myNgrains ! loop over grains
myInstance = phase_constitutionInstance(material_phase(g,i,e))
select case(phase_constitution(material_phase(g,i,e)))
@ -194,15 +195,15 @@ endif
allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
if (any(numerics_integrator == 1)) then
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
enddo
endif
@ -221,15 +222,15 @@ endif
allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
if (any(numerics_integrator == 1)) then
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
enddo
endif
@ -248,15 +249,15 @@ endif
allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
if (any(numerics_integrator == 1)) then
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
enddo
endif
@ -275,15 +276,15 @@ endif
allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
if (any(numerics_integrator == 1)) then
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
enddo
endif
@ -302,15 +303,15 @@ endif
allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
if (any(numerics_integrator == 1)) then
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 4)) then
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
endif
if (any(numerics_integrator == 5)) then
do s = 1,6
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
enddo
endif
@ -340,7 +341,7 @@ constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)
write(6,*) '<<<+- constitutive init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
if (debug_verbosity > 0) then
if (debug_verbosity > 0_pInt) then
write(6,'(a32,1x,7(i8,1x))') 'constitutive_state0: ', shape(constitutive_state0)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_subState0: ', shape(constitutive_subState0)
@ -457,11 +458,7 @@ endfunction
subroutine constitutive_microstructure(Temperature, Fe, Fp, ipc, ip, el)
use prec, only: pReal,pInt
use material, only: phase_constitution, &
material_phase, &
homogenization_maxNgrains
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_maxNipNeighbors
material_phase
use constitutive_j2, only: constitutive_j2_label, &
constitutive_j2_microstructure
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, &
@ -581,8 +578,7 @@ use debug, only: debug_cumDotStateCalls, &
debug_cumDotStateTicks, &
debug_verbosity
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_maxNipNeighbors
mesh_maxNips
use material, only: phase_constitution, &
material_phase, &
homogenization_maxNgrains
@ -620,7 +616,7 @@ integer(pLongInt) tick, tock, &
tickrate, &
maxticks
if (debug_verbosity > 0) then
if (debug_verbosity > 0_pInt) then
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif
@ -644,7 +640,7 @@ select case (phase_constitution(material_phase(ipc,ip,el)))
end select
if (debug_verbosity > 6) then
if (debug_verbosity > 6_pInt) then
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
!$OMP CRITICAL (debugTimingDotState)
debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt
@ -699,7 +695,7 @@ integer(pLongInt) tick, tock, &
maxticks
if (debug_verbosity > 0) then
if (debug_verbosity > 0_pInt) then
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif
@ -722,7 +718,7 @@ select case (phase_constitution(material_phase(ipc,ip,el)))
end select
if (debug_verbosity > 6) then
if (debug_verbosity > 6_pInt) then
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
!$OMP CRITICAL (debugTimingDotTemperature)
debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt
@ -748,8 +744,7 @@ function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el)
!*********************************************************************
use prec, only: pReal,pInt
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_maxNipNeighbors
mesh_maxNips
use material, only: phase_constitution, &
material_phase, &
homogenization_maxNgrains

398
code/constitutive_dislotwin.f90
View File

@ -133,7 +133,7 @@ subroutine constitutive_dislotwin_init(file)
!**************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt,pReal
use math, only: math_Mandel3333to66,math_Voigt66to3333,math_mul3x3
use mesh, only: mesh_maxNips, mesh_NcpElems
@ -144,9 +144,10 @@ use lattice
!* Input variables
integer(pInt), intent(in) :: file
!* Local variables
integer(pInt), parameter :: maxNchunks = 21
integer(pInt), parameter :: maxNchunks = 21_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section,maxNinstance,f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t1,t2,ns,nt,mySize,myStructure,maxTotalNslip,maxTotalNtwin
integer(pInt) :: section, maxNinstance,mySize,myStructure,maxTotalNslip,maxTotalNtwin,&
f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t1,t2,ns,nt
character(len=64) tag
character(len=1024) line
@ -157,130 +158,130 @@ character(len=1024) line
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
maxNinstance = count(phase_constitution == constitutive_dislotwin_label)
if (maxNinstance == 0) return
maxNinstance = int(count(phase_constitution == constitutive_dislotwin_label),pInt)
if (maxNinstance == 0_pInt) return
!* Space allocation for global variables
allocate(constitutive_dislotwin_sizeDotState(maxNinstance))
constitutive_dislotwin_sizeDotState = 0_pInt
allocate(constitutive_dislotwin_sizeState(maxNinstance))
constitutive_dislotwin_sizeState = 0_pInt
allocate(constitutive_dislotwin_sizePostResults(maxNinstance))
constitutive_dislotwin_sizePostResults = 0_pInt
allocate(constitutive_dislotwin_sizePostResult(maxval(phase_Noutput),maxNinstance))
constitutive_dislotwin_sizePostResult = 0_pInt
allocate(constitutive_dislotwin_output(maxval(phase_Noutput),maxNinstance))
constitutive_dislotwin_output = ''
allocate(constitutive_dislotwin_Noutput(maxNinstance))
constitutive_dislotwin_sizeDotState = 0_pInt
constitutive_dislotwin_sizeState = 0_pInt
constitutive_dislotwin_sizePostResults = 0_pInt
constitutive_dislotwin_sizePostResult = 0_pInt
constitutive_dislotwin_output = ''
constitutive_dislotwin_Noutput = 0_pInt
constitutive_dislotwin_Noutput = 0_pInt
allocate(constitutive_dislotwin_structureName(maxNinstance))
constitutive_dislotwin_structureName = ''
allocate(constitutive_dislotwin_structure(maxNinstance))
constitutive_dislotwin_structure = 0_pInt
allocate(constitutive_dislotwin_Nslip(lattice_maxNslipFamily,maxNinstance))
constitutive_dislotwin_Nslip = 0_pInt
allocate(constitutive_dislotwin_Ntwin(lattice_maxNtwinFamily,maxNinstance))
constitutive_dislotwin_Ntwin = 0_pInt
allocate(constitutive_dislotwin_slipFamily(lattice_maxNslip,maxNinstance))
constitutive_dislotwin_slipFamily = 0_pInt
allocate(constitutive_dislotwin_twinFamily(lattice_maxNtwin,maxNinstance))
constitutive_dislotwin_twinFamily = 0_pInt
allocate(constitutive_dislotwin_slipSystemLattice(lattice_maxNslip,maxNinstance))
constitutive_dislotwin_slipSystemLattice = 0_pInt
allocate(constitutive_dislotwin_twinSystemLattice(lattice_maxNtwin,maxNinstance))
constitutive_dislotwin_twinSystemLattice = 0_pInt
allocate(constitutive_dislotwin_totalNslip(maxNinstance))
constitutive_dislotwin_totalNslip = 0_pInt
allocate(constitutive_dislotwin_totalNtwin(maxNinstance))
constitutive_dislotwin_structureName = ''
constitutive_dislotwin_structure = 0_pInt
constitutive_dislotwin_Nslip = 0_pInt
constitutive_dislotwin_Ntwin = 0_pInt
constitutive_dislotwin_slipFamily = 0_pInt
constitutive_dislotwin_twinFamily = 0_pInt
constitutive_dislotwin_slipSystemLattice = 0_pInt
constitutive_dislotwin_twinSystemLattice = 0_pInt
constitutive_dislotwin_totalNslip = 0_pInt
constitutive_dislotwin_totalNtwin = 0_pInt
constitutive_dislotwin_totalNtwin = 0_pInt
allocate(constitutive_dislotwin_CoverA(maxNinstance))
constitutive_dislotwin_CoverA = 0.0_pReal
allocate(constitutive_dislotwin_C11(maxNinstance))
constitutive_dislotwin_C11 = 0.0_pReal
allocate(constitutive_dislotwin_C12(maxNinstance))
constitutive_dislotwin_C12 = 0.0_pReal
allocate(constitutive_dislotwin_C13(maxNinstance))
constitutive_dislotwin_C13 = 0.0_pReal
allocate(constitutive_dislotwin_C33(maxNinstance))
constitutive_dislotwin_C33 = 0.0_pReal
allocate(constitutive_dislotwin_C44(maxNinstance))
constitutive_dislotwin_C44 = 0.0_pReal
allocate(constitutive_dislotwin_Gmod(maxNinstance))
constitutive_dislotwin_Gmod = 0.0_pReal
allocate(constitutive_dislotwin_CAtomicVolume(maxNinstance))
constitutive_dislotwin_CAtomicVolume = 0.0_pReal
allocate(constitutive_dislotwin_D0(maxNinstance))
constitutive_dislotwin_D0 = 0.0_pReal
allocate(constitutive_dislotwin_Qsd(maxNinstance))
constitutive_dislotwin_Qsd = 0.0_pReal
allocate(constitutive_dislotwin_GrainSize(maxNinstance))
constitutive_dislotwin_GrainSize = 0.0_pReal
allocate(constitutive_dislotwin_p(maxNinstance))
constitutive_dislotwin_p = 0.0_pReal
allocate(constitutive_dislotwin_q(maxNinstance))
constitutive_dislotwin_q = 0.0_pReal
allocate(constitutive_dislotwin_MaxTwinFraction(maxNinstance))
constitutive_dislotwin_MaxTwinFraction = 0.0_pReal
allocate(constitutive_dislotwin_r(maxNinstance))
constitutive_dislotwin_r = 0.0_pReal
allocate(constitutive_dislotwin_CEdgeDipMinDistance(maxNinstance))
constitutive_dislotwin_CEdgeDipMinDistance = 0.0_pReal
allocate(constitutive_dislotwin_Cmfptwin(maxNinstance))
constitutive_dislotwin_Cmfptwin = 0.0_pReal
allocate(constitutive_dislotwin_Cthresholdtwin(maxNinstance))
constitutive_dislotwin_Cthresholdtwin = 0.0_pReal
allocate(constitutive_dislotwin_SolidSolutionStrength(maxNinstance))
constitutive_dislotwin_SolidSolutionStrength = 0.0_pReal
allocate(constitutive_dislotwin_L0(maxNinstance))
constitutive_dislotwin_L0 = 0.0_pReal
allocate(constitutive_dislotwin_aTolRho(maxNinstance))
constitutive_dislotwin_aTolRho = 0.0_pReal
allocate(constitutive_dislotwin_Cslip_66(6,6,maxNinstance))
constitutive_dislotwin_Cslip_66 = 0.0_pReal
allocate(constitutive_dislotwin_Cslip_3333(3,3,3,3,maxNinstance))
constitutive_dislotwin_Cslip_3333 = 0.0_pReal
allocate(constitutive_dislotwin_sbResistance(maxNinstance))
constitutive_dislotwin_sbResistance = 0.0_pReal
allocate(constitutive_dislotwin_sbVelocity(maxNinstance))
constitutive_dislotwin_sbVelocity = 0.0_pReal
allocate(constitutive_dislotwin_SFE_0K(maxNinstance))
constitutive_dislotwin_SFE_0K = 0.0_pReal
allocate(constitutive_dislotwin_dSFE_dT(maxNinstance))
constitutive_dislotwin_CoverA = 0.0_pReal
constitutive_dislotwin_C11 = 0.0_pReal
constitutive_dislotwin_C12 = 0.0_pReal
constitutive_dislotwin_C13 = 0.0_pReal
constitutive_dislotwin_C33 = 0.0_pReal
constitutive_dislotwin_C44 = 0.0_pReal
constitutive_dislotwin_Gmod = 0.0_pReal
constitutive_dislotwin_CAtomicVolume = 0.0_pReal
constitutive_dislotwin_D0 = 0.0_pReal
constitutive_dislotwin_Qsd = 0.0_pReal
constitutive_dislotwin_GrainSize = 0.0_pReal
constitutive_dislotwin_p = 0.0_pReal
constitutive_dislotwin_q = 0.0_pReal
constitutive_dislotwin_MaxTwinFraction = 0.0_pReal
constitutive_dislotwin_r = 0.0_pReal
constitutive_dislotwin_CEdgeDipMinDistance = 0.0_pReal
constitutive_dislotwin_Cmfptwin = 0.0_pReal
constitutive_dislotwin_Cthresholdtwin = 0.0_pReal
constitutive_dislotwin_SolidSolutionStrength= 0.0_pReal
constitutive_dislotwin_L0 = 0.0_pReal
constitutive_dislotwin_aTolRho = 0.0_pReal
constitutive_dislotwin_Cslip_66 = 0.0_pReal
constitutive_dislotwin_Cslip_3333 = 0.0_pReal
constitutive_dislotwin_sbResistance = 0.0_pReal
constitutive_dislotwin_sbVelocity = 0.0_pReal
constitutive_dislotwin_SFE_0K = 0.0_pReal
constitutive_dislotwin_dSFE_dT = 0.0_pReal
constitutive_dislotwin_dSFE_dT = 0.0_pReal
allocate(constitutive_dislotwin_rhoEdge0(lattice_maxNslipFamily,maxNinstance))
constitutive_dislotwin_rhoEdge0 = 0.0_pReal
allocate(constitutive_dislotwin_rhoEdgeDip0(lattice_maxNslipFamily,maxNinstance))
constitutive_dislotwin_rhoEdgeDip0 = 0.0_pReal
allocate(constitutive_dislotwin_burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_dislotwin_burgersPerSlipFamily = 0.0_pReal
allocate(constitutive_dislotwin_burgersPerTwinFamily(lattice_maxNtwinFamily,maxNinstance))
constitutive_dislotwin_burgersPerTwinFamily = 0.0_pReal
allocate(constitutive_dislotwin_QedgePerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_dislotwin_QedgePerSlipFamily = 0.0_pReal
allocate(constitutive_dislotwin_v0PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_dislotwin_v0PerSlipFamily = 0.0_pReal
allocate(constitutive_dislotwin_Ndot0PerTwinFamily(lattice_maxNtwinFamily,maxNinstance))
constitutive_dislotwin_Ndot0PerTwinFamily = 0.0_pReal
allocate(constitutive_dislotwin_twinsizePerTwinFamily(lattice_maxNtwinFamily,maxNinstance))
constitutive_dislotwin_twinsizePerTwinFamily = 0.0_pReal
allocate(constitutive_dislotwin_CLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_dislotwin_rhoEdge0 = 0.0_pReal
constitutive_dislotwin_rhoEdgeDip0 = 0.0_pReal
constitutive_dislotwin_burgersPerSlipFamily = 0.0_pReal
constitutive_dislotwin_burgersPerTwinFamily = 0.0_pReal
constitutive_dislotwin_QedgePerSlipFamily = 0.0_pReal
constitutive_dislotwin_v0PerSlipFamily = 0.0_pReal
constitutive_dislotwin_Ndot0PerTwinFamily = 0.0_pReal
constitutive_dislotwin_twinsizePerTwinFamily = 0.0_pReal
constitutive_dislotwin_CLambdaSlipPerSlipFamily = 0.0_pReal
constitutive_dislotwin_CLambdaSlipPerSlipFamily = 0.0_pReal
allocate(constitutive_dislotwin_interactionSlipSlip(lattice_maxNinteraction,maxNinstance))
constitutive_dislotwin_interactionSlipSlip = 0.0_pReal
allocate(constitutive_dislotwin_interactionSlipTwin(lattice_maxNinteraction,maxNinstance))
constitutive_dislotwin_interactionSlipTwin = 0.0_pReal
allocate(constitutive_dislotwin_interactionTwinSlip(lattice_maxNinteraction,maxNinstance))
constitutive_dislotwin_interactionTwinSlip = 0.0_pReal
allocate(constitutive_dislotwin_interactionTwinTwin(lattice_maxNinteraction,maxNinstance))
constitutive_dislotwin_interactionSlipSlip = 0.0_pReal
constitutive_dislotwin_interactionSlipTwin = 0.0_pReal
constitutive_dislotwin_interactionTwinSlip = 0.0_pReal
constitutive_dislotwin_interactionTwinTwin = 0.0_pReal
constitutive_dislotwin_interactionTwinTwin = 0.0_pReal
allocate(constitutive_dislotwin_sbSv(6,6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
constitutive_dislotwin_sbSv = 0.0_pReal
constitutive_dislotwin_sbSv = 0.0_pReal
!* Readout data from material.config file
rewind(file)
line = ''
section = 0
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
@ -297,7 +298,7 @@ do ! read thru sections of
if (section > 0_pInt .and. phase_constitution(section) == constitutive_dislotwin_label) then ! one of my sections
i = phase_constitutionInstance(section) ! which instance of my constitution is present phase
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('constitution')
cycle
@ -414,7 +415,7 @@ enddo
if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt,e=i)
if (sum(constitutive_dislotwin_Nslip(:,i)) <= 0_pInt) call IO_error(241_pInt,e=i,ext_msg='nslip')
if (sum(constitutive_dislotwin_Ntwin(:,i)) < 0_pInt) call IO_error(241_pInt,e=i,ext_msg='ntwin')
do f = 1,lattice_maxNslipFamily
do f = 1_pInt,lattice_maxNslipFamily
if (constitutive_dislotwin_Nslip(f,i) > 0_pInt) then
if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='rhoEdge0')
if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='rhoEdgeDip0')
@ -422,7 +423,7 @@ enddo
if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='v0')
endif
enddo
do f = 1,lattice_maxNtwinFamily
do f = 1_pInt,lattice_maxNtwinFamily
if (constitutive_dislotwin_Ntwin(f,i) > 0_pInt) then
if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='twinburgers')
if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='ndot0')
@ -450,35 +451,35 @@ maxTotalNslip = maxval(constitutive_dislotwin_totalNslip)
maxTotalNtwin = maxval(constitutive_dislotwin_totalNtwin)
allocate(constitutive_dislotwin_burgersPerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_dislotwin_burgersPerSlipSystem = 0.0_pReal
allocate(constitutive_dislotwin_burgersPerTwinSystem(maxTotalNtwin, maxNinstance))
constitutive_dislotwin_burgersPerTwinSystem= 0.0_pReal
allocate(constitutive_dislotwin_QedgePerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_dislotwin_QedgePerSlipSystem = 0.0_pReal
allocate(constitutive_dislotwin_v0PerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_dislotwin_v0PerSlipSystem = 0.0_pReal
allocate(constitutive_dislotwin_Ndot0PerTwinSystem(maxTotalNtwin, maxNinstance))
constitutive_dislotwin_Ndot0PerTwinSystem = 0.0_pReal
allocate(constitutive_dislotwin_twinsizePerTwinSystem(maxTotalNtwin, maxNinstance))
constitutive_dislotwin_twinsizePerTwinSystem = 0.0_pReal
allocate(constitutive_dislotwin_CLambdaSlipPerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_dislotwin_burgersPerSlipSystem = 0.0_pReal
constitutive_dislotwin_burgersPerTwinSystem = 0.0_pReal
constitutive_dislotwin_QedgePerSlipSystem = 0.0_pReal
constitutive_dislotwin_v0PerSlipSystem = 0.0_pReal
constitutive_dislotwin_Ndot0PerTwinSystem = 0.0_pReal
constitutive_dislotwin_twinsizePerTwinSystem = 0.0_pReal
constitutive_dislotwin_CLambdaSlipPerSlipSystem = 0.0_pReal
constitutive_dislotwin_CLambdaSlipPerSlipSystem = 0.0_pReal
allocate(constitutive_dislotwin_interactionMatrixSlipSlip(maxTotalNslip,maxTotalNslip,maxNinstance))
constitutive_dislotwin_interactionMatrixSlipSlip = 0.0_pReal
allocate(constitutive_dislotwin_interactionMatrixSlipTwin(maxTotalNslip,maxTotalNtwin,maxNinstance))
constitutive_dislotwin_interactionMatrixSlipTwin = 0.0_pReal
allocate(constitutive_dislotwin_interactionMatrixTwinSlip(maxTotalNtwin,maxTotalNslip,maxNinstance))
constitutive_dislotwin_interactionMatrixTwinSlip = 0.0_pReal
allocate(constitutive_dislotwin_interactionMatrixTwinTwin(maxTotalNtwin,maxTotalNtwin,maxNinstance))
constitutive_dislotwin_interactionMatrixTwinTwin = 0.0_pReal
allocate(constitutive_dislotwin_forestProjectionEdge(maxTotalNslip,maxTotalNslip,maxNinstance))
constitutive_dislotwin_interactionMatrixSlipSlip = 0.0_pReal
constitutive_dislotwin_interactionMatrixSlipTwin = 0.0_pReal
constitutive_dislotwin_interactionMatrixTwinSlip = 0.0_pReal
constitutive_dislotwin_interactionMatrixTwinTwin = 0.0_pReal
constitutive_dislotwin_forestProjectionEdge = 0.0_pReal
constitutive_dislotwin_forestProjectionEdge = 0.0_pReal
allocate(constitutive_dislotwin_Ctwin_66(6,6,maxTotalNtwin,maxNinstance))
constitutive_dislotwin_Ctwin_66 = 0.0_pReal
allocate(constitutive_dislotwin_Ctwin_3333(3,3,3,3,maxTotalNtwin,maxNinstance))
constitutive_dislotwin_Ctwin_66 = 0.0_pReal
constitutive_dislotwin_Ctwin_3333 = 0.0_pReal
constitutive_dislotwin_Ctwin_3333 = 0.0_pReal
do i = 1_pInt,maxNinstance
myStructure = constitutive_dislotwin_structure(i)
@ -504,11 +505,11 @@ do i = 1_pInt,maxNinstance
!* Determine size of state array
ns = constitutive_dislotwin_totalNslip(i)
nt = constitutive_dislotwin_totalNtwin(i)
constitutive_dislotwin_sizeDotState(i) = &
size(constitutive_dislotwin_listBasicSlipStates)*ns+size(constitutive_dislotwin_listBasicTwinStates)*nt
constitutive_dislotwin_sizeState(i) = &
constitutive_dislotwin_sizeDotState(i)+ &
size(constitutive_dislotwin_listDependentSlipStates)*ns+size(constitutive_dislotwin_listDependentTwinStates)*nt
constitutive_dislotwin_sizeDotState(i) =int(size(constitutive_dislotwin_listBasicSlipStates),pInt)*ns&
+int(size(constitutive_dislotwin_listBasicTwinStates),pInt)*nt
constitutive_dislotwin_sizeState(i) = constitutive_dislotwin_sizeDotState(i)&
+ int(size(constitutive_dislotwin_listDependentSlipStates),pInt)*ns&
+ int(size(constitutive_dislotwin_listDependentTwinStates),pInt)*nt
!* Determine size of postResults array
do o = 1_pInt,constitutive_dislotwin_Noutput(i)
@ -665,7 +666,7 @@ function constitutive_dislotwin_stateInit(myInstance)
!*********************************************************************
use prec, only: pReal,pInt
use math, only: pi
use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily
use lattice, only: lattice_maxNslipFamily
implicit none
!* Input-Output variables
@ -701,7 +702,7 @@ constitutive_dislotwin_stateInit(ns+1:2_pInt*ns) = rhoEdgeDip0
forall (s = 1_pInt:ns) &
invLambdaSlip0(s) = sqrt(dot_product((rhoEdge0+rhoEdgeDip0),constitutive_dislotwin_forestProjectionEdge(1:ns,s,myInstance)))/ &
constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,myInstance)
constitutive_dislotwin_stateInit(2_pInt*ns+nt+1:3_pInt*ns+nt) = invLambdaSlip0
constitutive_dislotwin_stateInit(2_pInt*ns+nt+1_pInt:3_pInt*ns+nt) = invLambdaSlip0
forall (s = 1_pInt:ns) &
MeanFreePathSlip0(s) = &
@ -715,14 +716,14 @@ sqrt(dot_product((rhoEdge0+rhoEdgeDip0),constitutive_dislotwin_interactionMatrix
constitutive_dislotwin_stateInit(5_pInt*ns+3_pInt*nt+1:6_pInt*ns+3_pInt*nt) = tauSlipThreshold0
!* Initialize dependent twin microstructural variables
forall (t = 1:nt) &
forall (t = 1_pInt:nt) &
MeanFreePathTwin0(t) = constitutive_dislotwin_GrainSize(myInstance)
constitutive_dislotwin_stateInit(5*ns+2*nt+1:5*ns+3*nt) = MeanFreePathTwin0
constitutive_dislotwin_stateInit(5_pInt*ns+2_pInt*nt+1_pInt:5_pInt*ns+3_pInt*nt) = MeanFreePathTwin0
forall (t = 1:nt) &
forall (t = 1_pInt:nt) &
TwinVolume0(t) = &
(pi/6.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)*MeanFreePathTwin0(t)**(2.0_pReal)
constitutive_dislotwin_stateInit(6*ns+4*nt+1:6*ns+5*nt) = TwinVolume0
constitutive_dislotwin_stateInit(6_pInt*ns+4_pInt*nt+1_pInt:6_pInt*ns+5_pInt*nt) = TwinVolume0
!write(6,*) '#STATEINIT#'
!write(6,*)
@ -782,13 +783,13 @@ ns = constitutive_dislotwin_totalNslip(myInstance)
nt = constitutive_dislotwin_totalNtwin(myInstance)
!* Total twin volume fraction
sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0
sumf = sum(state(g,ip,el)%p((2_pInt*ns+1_pInt):(2_pInt*ns+nt))) ! safe for nt == 0
!* Homogenized elasticity matrix
constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*constitutive_dislotwin_Cslip_66(:,:,myInstance)
do i=1,nt
do i=1_pInt,nt
constitutive_dislotwin_homogenizedC = &
constitutive_dislotwin_homogenizedC + state(g,ip,el)%p(2*ns+i)*constitutive_dislotwin_Ctwin_66(:,:,i,myInstance)
constitutive_dislotwin_homogenizedC + state(g,ip,el)%p(2_pInt*ns+i)*constitutive_dislotwin_Ctwin_66(:,:,i,myInstance)
enddo
return
@ -808,7 +809,6 @@ use prec, only: pReal,pInt,p_vec
use math, only: pi
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
use lattice, only: lattice_interactionSlipTwin,lattice_interactionTwinTwin
!use debug, only: debugger
implicit none
@ -846,69 +846,69 @@ sfe = constitutive_dislotwin_SFE_0K(myInstance) + &
constitutive_dislotwin_dSFE_dT(myInstance) * Temperature
!* rescaled twin volume fraction for topology
forall (t = 1:nt) &
forall (t = 1_pInt:nt) &
fOverStacksize(t) = &
state(g,ip,el)%p(2*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)
state(g,ip,el)%p(2_pInt*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)
!* 1/mean free distance between 2 forest dislocations seen by a moving dislocation
forall (s = 1:ns) &
state(g,ip,el)%p(2*ns+nt+s) = &
sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),&
forall (s = 1_pInt:ns) &
state(g,ip,el)%p(2_pInt*ns+nt+s) = &
sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1_pInt:2_pInt*ns)),&
constitutive_dislotwin_forestProjectionEdge(1:ns,s,myInstance)))/ &
constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,myInstance)
!* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
!$OMP CRITICAL (evilmatmul)
state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = 0.0_pReal
state(g,ip,el)%p((3_pInt*ns+nt+1_pInt):(4_pInt*ns+nt)) = 0.0_pReal
if (nt > 0_pInt) &
state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = &
state(g,ip,el)%p((3_pInt*ns+nt+1):(4_pInt*ns+nt)) = &
matmul(constitutive_dislotwin_interactionMatrixSlipTwin(1:ns,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf)
!$OMP END CRITICAL (evilmatmul)
!* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
!$OMP CRITICAL (evilmatmul)
if (nt > 0_pInt) &
state(g,ip,el)%p((4*ns+nt+1):(4*ns+2*nt)) = &
state(g,ip,el)%p((4_pInt*ns+nt+1_pInt):(4_pInt*ns+2_pInt*nt)) = &
matmul(constitutive_dislotwin_interactionMatrixTwinTwin(1:nt,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf)
!$OMP END CRITICAL (evilmatmul)
!* mean free path between 2 obstacles seen by a moving dislocation
do s = 1,ns
do s = 1_pInt,ns
if (nt > 0_pInt) then
state(g,ip,el)%p(4*ns+2*nt+s) = &
state(g,ip,el)%p(4_pInt*ns+2_pInt*nt+s) = &
constitutive_dislotwin_GrainSize(myInstance)/(1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*&
(state(g,ip,el)%p(2*ns+nt+s)+state(g,ip,el)%p(3*ns+nt+s)))
(state(g,ip,el)%p(2_pInt*ns+nt+s)+state(g,ip,el)%p(3_pInt*ns+nt+s)))
else
state(g,ip,el)%p(4*ns+s) = &
state(g,ip,el)%p(4_pInt*ns+s) = &
constitutive_dislotwin_GrainSize(myInstance)/&
(1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*(state(g,ip,el)%p(2*ns+s)))
(1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*(state(g,ip,el)%p(2_pInt*ns+s)))
endif
enddo
!* mean free path between 2 obstacles seen by a growing twin
forall (t = 1:nt) &
state(g,ip,el)%p(5*ns+2*nt+t) = &
forall (t = 1_pInt:nt) &
state(g,ip,el)%p(5_pInt*ns+2_pInt*nt+t) = &
(constitutive_dislotwin_Cmfptwin(myInstance)*constitutive_dislotwin_GrainSize(myInstance))/&
(1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*state(g,ip,el)%p(4*ns+nt+t))
(1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*state(g,ip,el)%p(4_pInt*ns+nt+t))
!* threshold stress for dislocation motion
forall (s = 1:ns) &
state(g,ip,el)%p(5*ns+3*nt+s) = constitutive_dislotwin_SolidSolutionStrength(myInstance)+ &
forall (s = 1_pInt:ns) &
state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+s) = constitutive_dislotwin_SolidSolutionStrength(myInstance)+ &
constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)*&
sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),&
sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1_pInt:2_pInt*ns)),&
constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance)))
!* threshold stress for growing twin
forall (t = 1:nt) &
state(g,ip,el)%p(6*ns+3*nt+t) = &
forall (t = 1_pInt:nt) &
state(g,ip,el)%p(6_pInt*ns+3_pInt*nt+t) = &
constitutive_dislotwin_Cthresholdtwin(myInstance)*&
(sfe/(3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,myInstance))+&
3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,myInstance)*constitutive_dislotwin_Gmod(myInstance)/&
(constitutive_dislotwin_L0(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(t,myInstance)))
!* final twin volume after growth
forall (t = 1:nt) &
state(g,ip,el)%p(6*ns+4*nt+t) = &
forall (t = 1_pInt:nt) &
state(g,ip,el)%p(6_pInt*ns+4_pInt*nt+t) = &
(pi/6.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)*state(g,ip,el)%p(5*ns+2*nt+t)**(2.0_pReal)
!if ((ip==1).and.(el==1)) then
@ -993,7 +993,7 @@ ns = constitutive_dislotwin_totalNslip(myInstance)
nt = constitutive_dislotwin_totalNtwin(myInstance)
!* Total twin volume fraction
sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0
sumf = sum(state(g,ip,el)%p((2_pInt*ns+1_pInt):(2_pInt*ns+nt))) ! safe for nt == 0
Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal
@ -1003,9 +1003,9 @@ dLp_dTstar = 0.0_pReal
gdot_slip = 0.0_pReal
dgdot_dtauslip = 0.0_pReal
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
j = j+1_pInt
!* Calculation of Lp
@ -1036,7 +1036,7 @@ do f = 1,lattice_maxNslipFamily ! loop over all
Lp = Lp + (1.0_pReal - sumf)*gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,myStructure)
!* Calculation of the tangent of Lp
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = &
dLp_dTstar3333(k,l,m,n) + dgdot_dtauslip(j)*&
lattice_Sslip(k,l,index_myFamily+i,myStructure)*&
@ -1049,7 +1049,7 @@ if(constitutive_dislotwin_sbVelocity(myInstance) /= 0.0_pReal) then
gdot_sb = 0.0_pReal
dgdot_dtausb = 0.0_pReal
call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
do j = 1,6
do j = 1_pInt,6_pInt
sb_s = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_sComposition(1:3,j))
sb_m = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_mComposition(1:3,j))
sb_Smatrix = math_tensorproduct(sb_s,sb_m)
@ -1073,20 +1073,20 @@ if(constitutive_dislotwin_sbVelocity(myInstance) /= 0.0_pReal) then
DotGamma0 = constitutive_dislotwin_sbVelocity(myInstance)
!* Shear rates due to shearband
gdot_sb(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*&
gdot_sb(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_q(myInstance))*&
sign(1.0_pReal,tau_sb(j))
!* Derivatives of shear rates
dgdot_dtausb(j) = &
((abs(gdot_sb(j))*BoltzmannRatio*&
constitutive_dislotwin_p(myInstance)*constitutive_dislotwin_q(myInstance))/constitutive_dislotwin_sbResistance(myInstance))*&
StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal)
StressRatio_pminus1*(1_pInt-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal)
!* Plastic velocity gradient for shear banding
Lp = Lp + gdot_sb(j)*sb_Smatrix
!* Calculation of the tangent of Lp
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = &
dLp_dTstar3333(k,l,m,n) + dgdot_dtausb(j)*&
sb_Smatrix(k,l)*&
@ -1098,9 +1098,9 @@ end if
gdot_twin = 0.0_pReal
dgdot_dtautwin = 0.0_pReal
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all slip families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all slip families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family
j = j+1_pInt
!* Calculation of Lp
@ -1122,7 +1122,7 @@ do f = 1,lattice_maxNtwinFamily ! loop over all
Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,myStructure)
!* Calculation of the tangent of Lp
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = &
dLp_dTstar3333(k,l,m,n) + dgdot_dtautwin(j)*&
lattice_Stwin(k,l,index_myFamily+i,myStructure)*&
@ -1164,11 +1164,9 @@ use prec, only: pReal,pInt,p_vec
use math, only: pi
use mesh, only: mesh_NcpElems, mesh_maxNips
use material, only: homogenization_maxNgrains, material_phase, phase_constitutionInstance
use lattice, only: lattice_Sslip,lattice_Sslip_v, &
lattice_Stwin,lattice_Stwin_v, &
use lattice, only: lattice_Sslip_v, lattice_Stwin_v, &
lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem, &
lattice_shearTwin
lattice_NslipSystem,lattice_NtwinSystem
implicit none
!* Input-Output variables
@ -1196,24 +1194,26 @@ ns = constitutive_dislotwin_totalNslip(myInstance)
nt = constitutive_dislotwin_totalNtwin(myInstance)
!* Total twin volume fraction
sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0
sumf = sum(state(g,ip,el)%p((2_pInt*ns+1_pInt):(2_pInt*ns+nt))) ! safe for nt == 0
constitutive_dislotwin_dotState = 0.0_pReal
!* Dislocation density evolution
gdot_slip = 0.0_pReal
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,MyStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,MyStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
j = j+1_pInt
!* Resolved shear stress on slip system
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
!* Stress ratios
StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance)
StressRatio_pminus1 = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal)
StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
constitutive_dislotwin_p(myInstance)
StressRatio_pminus1 = (abs(tau_slip(j))/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
(constitutive_dislotwin_p(myInstance)-1.0_pReal)
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature)
!* Initial shear rates
@ -1222,7 +1222,7 @@ do f = 1,lattice_maxNslipFamily ! loop over all
constitutive_dislotwin_v0PerSlipSystem(f,myInstance)
!* Shear rates due to slip
gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*&
gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_q(myInstance))*&
sign(1.0_pReal,tau_slip(j))
!* Multiplication
@ -1283,9 +1283,9 @@ enddo
!* Twin volume fraction evolution
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,MyStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) twin system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,MyStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) twin system in family
j = j+1_pInt
!* Resolved shear stress on twin system
@ -1295,9 +1295,9 @@ do f = 1,lattice_maxNtwinFamily ! loop over all
!* Shear rates and their derivatives due to twin
if ( tau_twin(j) > 0.0_pReal ) then
constitutive_dislotwin_dotState(2*ns+j) = &
constitutive_dislotwin_dotState(2_pInt*ns+j) = &
(constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*&
state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r)
state(g,ip,el)%p(6_pInt*ns+4_pInt*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r)
endif
enddo
@ -1365,7 +1365,7 @@ use math, only: pi,math_Mandel6to33, math_spectralDecompositionSym33
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput
use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin
lattice_NslipSystem,lattice_NtwinSystem
implicit none
!* Definition of variables
@ -1397,27 +1397,29 @@ constitutive_dislotwin_postResults = 0.0_pReal
!* Spectral decomposition of stress
call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
do o = 1,phase_Noutput(material_phase(g,ip,el))
do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
select case(constitutive_dislotwin_output(o,myInstance))
case ('edge_density')
constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns)
constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p(1:ns)
c = c + ns
case ('dipole_density')
constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p(ns+1:2*ns)
constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p(ns+1:2_pInt*ns)
c = c + ns
case ('shear_rate_slip')
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
j = j + 1_pInt
!* Resolved shear stress on slip system
tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
!* Stress ratios
StressRatio_p = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance)
StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal)
StressRatio_p = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
constitutive_dislotwin_p(myInstance)
StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
(constitutive_dislotwin_p(myInstance)-1.0_pReal)
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature)
!* Initial shear rates
@ -1427,29 +1429,33 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
!* Shear rates due to slip
constitutive_dislotwin_postResults(c+j) = &
DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau)
DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau)
enddo ; enddo
c = c + ns
case ('mfp_slip')
constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((4*ns+2*nt+1):(5*ns+2*nt))
constitutive_dislotwin_postResults(c+1_pInt:c+ns) =&
state(g,ip,el)%p((4_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+2_pInt*nt))
c = c + ns
case ('resolved_stress_slip')
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
j = j + 1_pInt
constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
constitutive_dislotwin_postResults(c+j) =&
dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
enddo; enddo
c = c + ns
case ('threshold_stress_slip')
constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((5*ns+3*nt+1):(6*ns+3*nt))
constitutive_dislotwin_postResults(c+1_pInt:c+ns) = &
state(g,ip,el)%p((5_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+3_pInt*nt))
c = c + ns
case ('edge_dipole_distance')
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
j = j + 1_pInt
constitutive_dislotwin_postResults(c+j) = &
(3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/&
@ -1459,15 +1465,15 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
enddo; enddo
c = c + ns
case ('resolved_stress_shearband')
do j = 1,6 ! loop over all shearband families
do j = 1_pInt,6_pInt ! loop over all shearband families
constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v, constitutive_dislotwin_sbSv(1:6,j,g,ip,el))
enddo
c = c + 6
c = c + 6_pInt
case ('schmid_factor_shearband')
constitutive_dislotwin_postResults(c+1:c+6) = constitutive_dislotwin_sbSv(1:6,j,g,ip,el)
c = c + 6
constitutive_dislotwin_postResults(c+1_pInt:c+6_pInt) = constitutive_dislotwin_sbSv(1:6,j,g,ip,el)
c = c + 6_pInt
case ('shear_rate_shearband')
do j = 1,6 ! loop over all shearband families
do j = 1_pInt,6_pInt ! loop over all shearband families
!* Resolved shear stress on shearband system
tau = dot_product(Tstar_v,constitutive_dislotwin_sbSv(1:6,j,g,ip,el))
!* Stress ratios
@ -1481,64 +1487,66 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
!* Shear rates due to slip
constitutive_dislotwin_postResults(c+j) = &
DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau)
DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau)
enddo
c = c + 6
c = c + 6_pInt
case ('twin_fraction')
constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((2*ns+1):(2*ns+nt))
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(g,ip,el)%p((2_pInt*ns+1_pInt):(2_pInt*ns+nt))
c = c + nt
case ('shear_rate_twin')
if (nt > 0_pInt) then
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) twin system in family
j = j + 1_pInt
!* Resolved shear stress on twin system
tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure))
!* Stress ratios
StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau)**constitutive_dislotwin_r(myInstance)
StressRatio_r = (state(g,ip,el)%p(6_pInt*ns+3_pInt*nt+j)/tau)**constitutive_dislotwin_r(myInstance)
!* Shear rates and their derivatives due to twin
if ( tau > 0.0_pReal ) then
constitutive_dislotwin_postResults(c+j) = &
(constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*&
state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r)
state(g,ip,el)%p(6_pInt*ns+4_pInt*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r)
endif
enddo ; enddo
endif
c = c + nt
case ('mfp_twin')
constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((5*ns+2*nt+1):(5*ns+3*nt))
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(g,ip,el)%p((5_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+3_pInt*nt))
c = c + nt
case ('resolved_stress_twin')
if (nt > 0_pInt) then
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all slip families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all slip families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family
j = j + 1_pInt
constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure))
enddo; enddo
endif
c = c + nt
case ('threshold_stress_twin')
constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((6*ns+3*nt+1):(6*ns+4*nt))
constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(g,ip,el)%p((6_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+4_pInt*nt))
c = c + nt
case ('stress_exponent')
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
j = j + 1_pInt
!* Resolved shear stress on slip system
tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
!* Stress ratios
StressRatio_p = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance)
StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal)
StressRatio_p = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
constitutive_dislotwin_p(myInstance)
StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
(constitutive_dislotwin_p(myInstance)-1.0_pReal)
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature)
!* Initial shear rates
@ -1547,14 +1555,14 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
constitutive_dislotwin_v0PerSlipSystem(f,myInstance)
!* Shear rates due to slip
gdot_slip = &
DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau)
gdot_slip = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau)
!* Derivatives of shear rates
dgdot_dtauslip = &
((abs(gdot_slip)*BoltzmannRatio*&
constitutive_dislotwin_p(myInstance)*constitutive_dislotwin_q(myInstance))/state(g,ip,el)%p(5*ns+3*nt+j))*&
StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal)
StressRatio_pminus1*(1_pInt-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal)
!* Stress exponent
if (gdot_slip==0.0_pReal) then
@ -1565,12 +1573,12 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
enddo ; enddo
c = c + ns
case ('sb_eigenvalues')
forall (j = 1:3) &
forall (j = 1_pInt:3_pInt) &
constitutive_dislotwin_postResults(c+j) = eigValues(j)
c = c + 3
c = c + 3_pInt
case ('sb_eigenvectors')
constitutive_dislotwin_postResults(c+1:c+9) = reshape(eigVectors,(/9/))
c = c + 9
constitutive_dislotwin_postResults(c+1_pInt:c+9_pInt) = reshape(eigVectors,(/9/))
c = c + 9_pInt
end select
enddo

77
code/constitutive_j2.f90
View File

@ -84,14 +84,15 @@ subroutine constitutive_j2_init(file)
!**************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333
use IO
use material
use debug, only: debug_verbosity
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt), parameter :: maxNchunks = 7_pInt
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k, mySize
character(len=64) tag
character(len=1024) line
@ -103,10 +104,10 @@ subroutine constitutive_j2_init(file)
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
maxNinstance = count(phase_constitution == constitutive_j2_label)
if (maxNinstance == 0) return
maxNinstance = int(count(phase_constitution == constitutive_j2_label),pInt)
if (maxNinstance == 0_pInt) return
if (debug_verbosity > 0) then
if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a16,1x,i5)') '# instances:',maxNinstance
write(6,*)
@ -133,7 +134,7 @@ subroutine constitutive_j2_init(file)
rewind(file)
line = ''
section = 0
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
@ -156,34 +157,34 @@ subroutine constitutive_j2_init(file)
cycle
case ('(output)')
constitutive_j2_Noutput(i) = constitutive_j2_Noutput(i) + 1_pInt
constitutive_j2_output(constitutive_j2_Noutput(i),i) = IO_lc(IO_stringValue(line,positions,2))
constitutive_j2_output(constitutive_j2_Noutput(i),i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('c11')
constitutive_j2_C11(i) = IO_floatValue(line,positions,2)
constitutive_j2_C11(i) = IO_floatValue(line,positions,2_pInt)
case ('c12')
constitutive_j2_C12(i) = IO_floatValue(line,positions,2)
constitutive_j2_C12(i) = IO_floatValue(line,positions,2_pInt)
case ('tau0')
constitutive_j2_tau0(i) = IO_floatValue(line,positions,2)
constitutive_j2_tau0(i) = IO_floatValue(line,positions,2_pInt)
case ('gdot0')
constitutive_j2_gdot0(i) = IO_floatValue(line,positions,2)
constitutive_j2_gdot0(i) = IO_floatValue(line,positions,2_pInt)
case ('n')
constitutive_j2_n(i) = IO_floatValue(line,positions,2)
constitutive_j2_n(i) = IO_floatValue(line,positions,2_pInt)
case ('h0')
constitutive_j2_h0(i) = IO_floatValue(line,positions,2)
constitutive_j2_h0(i) = IO_floatValue(line,positions,2_pInt)
case ('tausat')
constitutive_j2_tausat(i) = IO_floatValue(line,positions,2)
constitutive_j2_tausat(i) = IO_floatValue(line,positions,2_pInt)
case ('a', 'w0')
constitutive_j2_a(i) = IO_floatValue(line,positions,2)
constitutive_j2_a(i) = IO_floatValue(line,positions,2_pInt)
case ('taylorfactor')
constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2)
constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2_pInt)
case ('atol_resistance')
constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2)
constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2_pInt)
case default
call IO_error(210_pInt,ext_msg=tag)
end select
endif
enddo
100 do i = 1,maxNinstance ! sanity checks
100 do i = 1_pInt,maxNinstance ! sanity checks
if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(211_pInt,ext_msg='tau0')
if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='gdot0')
if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='n')
@ -193,8 +194,8 @@ subroutine constitutive_j2_init(file)
if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='aTol_resistance')
enddo
do i = 1,maxNinstance
do j = 1,constitutive_j2_Noutput(i)
do i = 1_pInt,maxNinstance
do j = 1_pInt,constitutive_j2_Noutput(i)
select case(constitutive_j2_output(j,i))
case('flowstress')
mySize = 1_pInt
@ -211,13 +212,14 @@ subroutine constitutive_j2_init(file)
endif
enddo
constitutive_j2_sizeDotState(i) = 1
constitutive_j2_sizeState(i) = 1
constitutive_j2_sizeDotState(i) = 1_pInt
constitutive_j2_sizeState(i) = 1_pInt
forall(k=1:3)
forall(j=1:3) &
constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i)
constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i)
forall(k=1_pInt:3_pInt)
forall(j=1_pInt:3_pInt)
constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i)
end forall
constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i)
constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i))
end forall
constitutive_j2_Cslip_66(1:6,1:6,i) = &
@ -333,8 +335,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v,
use math, only: math_mul6x6, &
math_Mandel6to33, &
math_Plain3333to99
use lattice, only: lattice_Sslip, &
lattice_Sslip_v
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
@ -379,7 +379,7 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v,
dLp_dTstar_99 = 0.0_pReal
! for Tstar==0 both Lp and dLp_dTstar are zero (if not n==1)
if (norm_Tstar_dev > 0) then
if (norm_Tstar_dev > 0_pInt) then
! Calculation of gamma_dot
gamma_dot = constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
@ -390,9 +390,9 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v,
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/constitutive_j2_fTaylor(matID)
!* Calculation of the tangent of Lp
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,m,n) = (constitutive_j2_n(matID)-1.0_pReal) * Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
forall (k=1:3,l=1:3) &
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
dLp_dTstar_99 = math_Plain3333to99(gamma_dot / constitutive_j2_fTaylor(matID) * dLp_dTstar_3333 / norm_Tstar_dev)
end if
@ -412,7 +412,6 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el)
pInt, &
p_vec
use math, only: math_mul6x6
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
@ -471,9 +470,8 @@ pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip,
!*** variables and functions from other modules ***!
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
use material, only: homogenization_maxNgrains
implicit none
!*** input variables ***!
@ -504,7 +502,6 @@ pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip
pInt, &
p_vec
use math, only: math_mul6x6
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
@ -550,17 +547,17 @@ pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip
c = 0_pInt
constitutive_j2_postResults = 0.0_pReal
do o = 1,phase_Noutput(material_phase(g,ip,el))
do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
select case(constitutive_j2_output(o,matID))
case ('flowstress')
constitutive_j2_postResults(c+1) = state(g,ip,el)%p(1)
c = c + 1
constitutive_j2_postResults(c+1_pInt) = state(g,ip,el)%p(1)
c = c + 1_pInt
case ('strainrate')
constitutive_j2_postResults(c+1) = &
constitutive_j2_postResults(c+1_pInt) = &
constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!---------------------------------------------------
(constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) ** constitutive_j2_n(matID)
c = c + 1
c = c + 1_pInt
end select
enddo

114
code/constitutive_phenopowerlaw.f90
View File

@ -145,10 +145,11 @@ CONTAINS
!****************************************
subroutine constitutive_phenopowerlaw_init(file)
subroutine constitutive_phenopowerlaw_init(myFile)
!**************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333
use IO
@ -163,7 +164,7 @@ subroutine constitutive_phenopowerlaw_init(file)
lattice_interactionTwinSlip, &
lattice_interactionTwinTwin
integer(pInt), intent(in) :: file
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = lattice_maxNinteraction + 1_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k, f,o, &
@ -178,7 +179,7 @@ subroutine constitutive_phenopowerlaw_init(file)
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
maxNinstance = count(phase_constitution == constitutive_phenopowerlaw_label)
maxNinstance = int(count(phase_constitution == constitutive_phenopowerlaw_label),pInt)
if (maxNinstance == 0) return
if (debug_verbosity > 0) then
@ -252,16 +253,16 @@ subroutine constitutive_phenopowerlaw_init(file)
allocate(constitutive_phenopowerlaw_aTolResistance(maxNinstance))
constitutive_phenopowerlaw_aTolResistance = 0.0_pReal
rewind(file)
rewind(myFile)
line = ''
section = 0
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
read(myFile,'(a1024)',END=100) line
enddo
do ! read thru sections of phase part
read(file,'(a1024)',END=100) line
read(myFile,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
@ -542,14 +543,14 @@ function constitutive_phenopowerlaw_stateInit(myInstance)
constitutive_phenopowerlaw_stateInit = 0.0_pReal
do i = 1,lattice_maxNslipFamily
do i = 1_pInt,lattice_maxNslipFamily
constitutive_phenopowerlaw_stateInit(1+&
sum(constitutive_phenopowerlaw_Nslip(1:i-1,myInstance)) : &
sum(constitutive_phenopowerlaw_Nslip(1:i ,myInstance))) = &
constitutive_phenopowerlaw_tau0_slip(i,myInstance)
enddo
do i = 1,lattice_maxNtwinFamily
do i = 1_pInt,lattice_maxNtwinFamily
constitutive_phenopowerlaw_stateInit(1+sum(constitutive_phenopowerlaw_Nslip(:,myInstance))+&
sum(constitutive_phenopowerlaw_Ntwin(1:i-1,myInstance)) : &
sum(constitutive_phenopowerlaw_Nslip(:,myInstance))+&
@ -677,17 +678,17 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
nSlip = constitutive_phenopowerlaw_totalNslip(matID)
nTwin = constitutive_phenopowerlaw_totalNtwin(matID)
index_Gamma = nSlip + nTwin + 1
index_F = nSlip + nTwin + 2
index_Gamma = nSlip + nTwin + 1_pInt
index_F = nSlip + nTwin + 2_pInt
Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal
dLp_dTstar = 0.0_pReal
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt
!* Calculation of Lp
@ -702,7 +703,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
if (gdot_slip(j) /= 0.0_pReal) then
dgdot_dtauslip(j) = gdot_slip(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip(j)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip(j)*lattice_Sslip(k,l,index_myFamily+i,structID)* &
lattice_Sslip(m,n,index_myFamily+i,structID)
@ -711,9 +712,9 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
enddo
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt
!* Calculation of Lp
@ -729,7 +730,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
if (gdot_twin(j) /= 0.0_pReal) then
dgdot_dtautwin(j) = gdot_twin(j)*constitutive_phenopowerlaw_n_twin(matID)/tau_twin(j)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtautwin(j)*lattice_Stwin(k,l,index_myFamily+i,structID)* &
lattice_Stwin(m,n,index_myFamily+i,structID)
@ -755,7 +756,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
!* - constitutive_dotState : evolution of state variable *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
use lattice, only: lattice_Sslip_v, lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
@ -780,8 +781,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
nSlip = constitutive_phenopowerlaw_totalNslip(matID)
nTwin = constitutive_phenopowerlaw_totalNtwin(matID)
index_Gamma = nSlip + nTwin + 1
index_F = nSlip + nTwin + 2
index_Gamma = nSlip + nTwin + 1_pInt
index_F = nSlip + nTwin + 2_pInt
constitutive_phenopowerlaw_dotState = 0.0_pReal
@ -800,9 +801,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
ssat_offset = constitutive_phenopowerlaw_spr(matID)*sqrt(state(ipc,ip,el)%p(index_F))
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt
h_slipslip(j) = c_slipslip*(1.0_pReal-state(ipc,ip,el)%p(j) / & ! system-dependent prefactor for slip--slip interaction
(constitutive_phenopowerlaw_tausat_slip(f,matID)+ssat_offset))** &
@ -819,9 +820,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
enddo
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt
h_twinslip(j) = c_twinslip ! no system-dependent parts
h_twintwin(j) = c_twintwin
@ -839,8 +840,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
!-- calculate the overall hardening based on above
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt
constitutive_phenopowerlaw_dotState(j) = & ! evolution of slip resistance j
h_slipslip(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_slipslip(1:nSlip,j,matID),abs(gdot_slip)) + & ! dot gamma_slip
@ -851,9 +852,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
enddo
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt
constitutive_phenopowerlaw_dotState(j+nSlip) = & ! evolution of twin resistance j
h_twinslip(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_twinslip(1:nSlip,j,matID),abs(gdot_slip)) + & ! dot gamma_slip
@ -875,9 +876,8 @@ pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,stat
!*** variables and functions from other modules ***!
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
use mesh, only: mesh_NcpElems, mesh_maxNips
use material, only: homogenization_maxNgrains
implicit none
!*** input variables ***!
@ -932,23 +932,23 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
nSlip = constitutive_phenopowerlaw_totalNslip(matID)
nTwin = constitutive_phenopowerlaw_totalNtwin(matID)
index_Gamma = nSlip + nTwin + 1
index_F = nSlip + nTwin + 2
index_Gamma = nSlip + nTwin + 1_pInt
index_F = nSlip + nTwin + 2_pInt
constitutive_phenopowerlaw_postResults = 0.0_pReal
c = 0_pInt
do o = 1,phase_Noutput(material_phase(ipc,ip,el))
do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_phenopowerlaw_output(o,matID))
case ('resistance_slip')
constitutive_phenopowerlaw_postResults(c+1:c+nSlip) = state(ipc,ip,el)%p(1:nSlip)
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(ipc,ip,el)%p(1:nSlip)
c = c + nSlip
case ('shearrate_slip')
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j + 1_pInt
tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID))
constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(matID)*&
@ -959,27 +959,27 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
case ('resolvedstress_slip')
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID))
enddo; enddo
c = c + nSlip
case ('totalshear')
constitutive_phenopowerlaw_postResults(c+1) = state(ipc,ip,el)%p(index_Gamma)
c = c + 1
constitutive_phenopowerlaw_postResults(c+1_pInt) = state(ipc,ip,el)%p(index_Gamma)
c = c + 1_pInt
case ('resistance_twin')
constitutive_phenopowerlaw_postResults(c+1:c+nTwin) = state(ipc,ip,el)%p(1+nSlip:nTwin+nSlip)
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = state(ipc,ip,el)%p(1_pInt+nSlip:nTwin+nSlip)
c = c + nTwin
case ('shearrate_twin')
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j + 1_pInt
tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
@ -991,17 +991,17 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
case ('resolvedstress_twin')
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
enddo; enddo
c = c + nTwin
case ('totalvolfrac')
constitutive_phenopowerlaw_postResults(c+1) = state(ipc,ip,el)%p(index_F)
c = c + 1
constitutive_phenopowerlaw_postResults(c+1_pInt) = state(ipc,ip,el)%p(index_F)
c = c + 1_pInt
end select
enddo

443
code/constitutive_titanmod.f90
View File

@ -203,7 +203,7 @@ subroutine constitutive_titanmod_init(file)
!**************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt,pReal
use math, only: math_Mandel3333to66,math_Voigt66to3333,math_mul3x3
use IO
@ -213,10 +213,11 @@ use lattice
!* Input variables
integer(pInt), intent(in) :: file
!* Local variables
integer(pInt), parameter :: maxNchunks = 21
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section,maxNinstance,f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t,t1,t2,ns,nt,mySize,myStructure,maxTotalNslip, &
integer(pInt), parameter :: maxNchunks = 21_pInt
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
integer(pInt) section,f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t,t1,t2,ns,nt,mySize,myStructure,maxTotalNslip, &
maxTotalNtwin
integer :: maxNinstance !no pInt
character(len=64) tag
character(len=1024) line
@ -229,158 +230,156 @@ maxNinstance = count(phase_constitution == constitutive_titanmod_label)
if (maxNinstance == 0) return
!* Space allocation for global variables
allocate(constitutive_titanmod_sizeDotState(maxNinstance))
allocate(constitutive_titanmod_sizeDotState(maxNinstance))
constitutive_titanmod_sizeDotState = 0_pInt
allocate(constitutive_titanmod_sizeState(maxNinstance))
constitutive_titanmod_sizeState = 0_pInt
allocate(constitutive_titanmod_sizePostResults(maxNinstance))
constitutive_titanmod_sizePostResults = 0_pInt
allocate(constitutive_titanmod_sizePostResult(maxval(phase_Noutput),maxNinstance))
constitutive_titanmod_sizePostResult = 0_pInt
allocate(constitutive_titanmod_output(maxval(phase_Noutput),maxNinstance))
constitutive_titanmod_output = ''
allocate(constitutive_titanmod_Noutput(maxNinstance))
constitutive_titanmod_sizeDotState = 0_pInt
constitutive_titanmod_sizeState = 0_pInt
constitutive_titanmod_sizePostResults = 0_pInt
constitutive_titanmod_sizePostResult = 0_pInt
constitutive_titanmod_output = ''
constitutive_titanmod_Noutput = 0_pInt
constitutive_titanmod_Noutput = 0_pInt
allocate(constitutive_titanmod_structureName(maxNinstance))
constitutive_titanmod_structureName = ''
allocate(constitutive_titanmod_structure(maxNinstance))
constitutive_titanmod_structure = 0_pInt
allocate(constitutive_titanmod_Nslip(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_Nslip = 0_pInt
allocate(constitutive_titanmod_Ntwin(lattice_maxNtwinFamily,maxNinstance))
constitutive_titanmod_Ntwin = 0_pInt
allocate(constitutive_titanmod_slipFamily(lattice_maxNslip,maxNinstance))
constitutive_titanmod_slipFamily = 0_pInt
allocate(constitutive_titanmod_twinFamily(lattice_maxNtwin,maxNinstance))
constitutive_titanmod_twinFamily = 0_pInt
allocate(constitutive_titanmod_slipSystemLattice(lattice_maxNslip,maxNinstance))
constitutive_titanmod_slipSystemLattice = 0_pInt
allocate(constitutive_titanmod_twinSystemLattice(lattice_maxNtwin,maxNinstance))
constitutive_titanmod_twinSystemLattice = 0_pInt
allocate(constitutive_titanmod_totalNslip(maxNinstance))
constitutive_titanmod_totalNslip = 0_pInt
allocate(constitutive_titanmod_totalNtwin(maxNinstance))
constitutive_titanmod_structureName = ''
constitutive_titanmod_structure = 0_pInt
constitutive_titanmod_Nslip = 0_pInt
constitutive_titanmod_Ntwin = 0_pInt
constitutive_titanmod_slipFamily = 0_pInt
constitutive_titanmod_twinFamily = 0_pInt
constitutive_titanmod_slipSystemLattice = 0_pInt
constitutive_titanmod_twinSystemLattice = 0_pInt
constitutive_titanmod_totalNslip = 0_pInt
constitutive_titanmod_totalNtwin = 0_pInt
constitutive_titanmod_totalNtwin = 0_pInt
allocate(constitutive_titanmod_CoverA(maxNinstance))
constitutive_titanmod_CoverA = 0.0_pReal
allocate(constitutive_titanmod_C11(maxNinstance))
constitutive_titanmod_C11 = 0.0_pReal
allocate(constitutive_titanmod_C12(maxNinstance))
constitutive_titanmod_C12 = 0.0_pReal
allocate(constitutive_titanmod_C13(maxNinstance))
constitutive_titanmod_C13 = 0.0_pReal
allocate(constitutive_titanmod_C33(maxNinstance))
constitutive_titanmod_C33 = 0.0_pReal
allocate(constitutive_titanmod_C44(maxNinstance))
constitutive_titanmod_C44 = 0.0_pReal
allocate(constitutive_titanmod_debyefrequency(maxNinstance))
constitutive_titanmod_debyefrequency = 0.0_pReal
allocate(constitutive_titanmod_kinkf0(maxNinstance))
constitutive_titanmod_kinkf0 = 0.0_pReal
allocate(constitutive_titanmod_Gmod(maxNinstance))
constitutive_titanmod_Gmod = 0.0_pReal
allocate(constitutive_titanmod_CAtomicVolume(maxNinstance))
constitutive_titanmod_CAtomicVolume = 0.0_pReal
allocate(constitutive_titanmod_dc(maxNinstance))
constitutive_titanmod_dc = 0.0_pReal
allocate(constitutive_titanmod_twinhpconstant(maxNinstance))
constitutive_titanmod_twinhpconstant = 0.0_pReal
allocate(constitutive_titanmod_GrainSize(maxNinstance))
constitutive_titanmod_GrainSize = 0.0_pReal
allocate(constitutive_titanmod_MaxTwinFraction(maxNinstance))
constitutive_titanmod_MaxTwinFraction = 0.0_pReal
allocate(constitutive_titanmod_r(maxNinstance))
constitutive_titanmod_r = 0.0_pReal
allocate(constitutive_titanmod_CEdgeDipMinDistance(maxNinstance))
constitutive_titanmod_CEdgeDipMinDistance = 0.0_pReal
allocate(constitutive_titanmod_Cmfptwin(maxNinstance))
constitutive_titanmod_Cmfptwin = 0.0_pReal
allocate(constitutive_titanmod_Cthresholdtwin(maxNinstance))
constitutive_titanmod_Cthresholdtwin = 0.0_pReal
allocate(constitutive_titanmod_aTolRho(maxNinstance))
constitutive_titanmod_aTolRho = 0.0_pReal
allocate(constitutive_titanmod_Cslip_66(6,6,maxNinstance))
constitutive_titanmod_Cslip_66 = 0.0_pReal
allocate(constitutive_titanmod_Cslip_3333(3,3,3,3,maxNinstance))
constitutive_titanmod_CoverA = 0.0_pReal
constitutive_titanmod_C11 = 0.0_pReal
constitutive_titanmod_C12 = 0.0_pReal
constitutive_titanmod_C13 = 0.0_pReal
constitutive_titanmod_C33 = 0.0_pReal
constitutive_titanmod_C44 = 0.0_pReal
constitutive_titanmod_debyefrequency = 0.0_pReal
constitutive_titanmod_kinkf0 = 0.0_pReal
constitutive_titanmod_Gmod = 0.0_pReal
constitutive_titanmod_CAtomicVolume = 0.0_pReal
constitutive_titanmod_dc = 0.0_pReal
constitutive_titanmod_twinhpconstant = 0.0_pReal
constitutive_titanmod_GrainSize = 0.0_pReal
constitutive_titanmod_MaxTwinFraction = 0.0_pReal
constitutive_titanmod_r = 0.0_pReal
constitutive_titanmod_CEdgeDipMinDistance = 0.0_pReal
constitutive_titanmod_Cmfptwin = 0.0_pReal
constitutive_titanmod_Cthresholdtwin = 0.0_pReal
constitutive_titanmod_aTolRho = 0.0_pReal
constitutive_titanmod_Cslip_66 = 0.0_pReal
constitutive_titanmod_Cslip_3333 = 0.0_pReal
constitutive_titanmod_Cslip_3333 = 0.0_pReal
allocate(constitutive_titanmod_rho_edge0(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_rho_edge0 = 0.0_pReal
allocate(constitutive_titanmod_rho_screw0(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_rho_screw0 = 0.0_pReal
allocate(constitutive_titanmod_shear_system0(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_shear_system0 = 0.0_pReal
allocate(constitutive_titanmod_burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_burgersPerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_burgersPerTwinFamily(lattice_maxNtwinFamily,maxNinstance))
constitutive_titanmod_burgersPerTwinFamily = 0.0_pReal
allocate(constitutive_titanmod_f0_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_f0_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_tau0e_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_tau0e_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_tau0s_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_tau0s_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_capre_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_capre_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_caprs_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_caprs_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_pe_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_pe_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_ps_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_ps_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_qe_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_qe_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_qs_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_qs_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_v0e_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_v0e_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_v0s_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_v0s_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_kinkcriticallength_PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_kinkcriticallength_PerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_twinsizePerTwinFamily(lattice_maxNtwinFamily,maxNinstance))
constitutive_titanmod_twinsizePerTwinFamily = 0.0_pReal
allocate(constitutive_titanmod_CeLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_CeLambdaSlipPerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_CsLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance))
constitutive_titanmod_CsLambdaSlipPerSlipFamily = 0.0_pReal
allocate(constitutive_titanmod_twinf0_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance))
constitutive_titanmod_twinf0_PerTwinFamily = 0.0_pReal
allocate(constitutive_titanmod_twinshearconstant_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance))
constitutive_titanmod_twinshearconstant_PerTwinFamily = 0.0_pReal
allocate(constitutive_titanmod_twintau0_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance))
constitutive_titanmod_twintau0_PerTwinFamily = 0.0_pReal
allocate(constitutive_titanmod_twinp_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance))
constitutive_titanmod_twingamma0_PerTwinFamily = 0.0_pReal
allocate(constitutive_titanmod_twinq_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance))
constitutive_titanmod_twinLambdaSlipPerTwinFamily = 0.0_pReal
allocate(constitutive_titanmod_twingamma0_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance))
constitutive_titanmod_twinp_PerTwinFamily = 0.0_pReal
allocate(constitutive_titanmod_twinLambdaSlipPerTwinFamily(lattice_maxNTwinFamily,maxNinstance))
constitutive_titanmod_rho_edge0 = 0.0_pReal
constitutive_titanmod_rho_screw0 = 0.0_pReal
constitutive_titanmod_shear_system0 = 0.0_pReal
constitutive_titanmod_burgersPerSlipFamily = 0.0_pReal
constitutive_titanmod_burgersPerTwinFamily = 0.0_pReal
constitutive_titanmod_f0_PerSlipFamily = 0.0_pReal
constitutive_titanmod_tau0e_PerSlipFamily = 0.0_pReal
constitutive_titanmod_tau0s_PerSlipFamily = 0.0_pReal
constitutive_titanmod_capre_PerSlipFamily = 0.0_pReal
constitutive_titanmod_caprs_PerSlipFamily = 0.0_pReal
constitutive_titanmod_v0e_PerSlipFamily = 0.0_pReal
constitutive_titanmod_v0s_PerSlipFamily = 0.0_pReal
constitutive_titanmod_kinkcriticallength_PerSlipFamily = 0.0_pReal
constitutive_titanmod_twinsizePerTwinFamily = 0.0_pReal
constitutive_titanmod_CeLambdaSlipPerSlipFamily = 0.0_pReal
constitutive_titanmod_CsLambdaSlipPerSlipFamily = 0.0_pReal
constitutive_titanmod_pe_PerSlipFamily = 0.0_pReal
constitutive_titanmod_ps_PerSlipFamily = 0.0_pReal
constitutive_titanmod_qe_PerSlipFamily = 0.0_pReal
constitutive_titanmod_qs_PerSlipFamily = 0.0_pReal
constitutive_titanmod_twinf0_PerTwinFamily = 0.0_pReal
constitutive_titanmod_twinshearconstant_PerTwinFamily = 0.0_pReal
constitutive_titanmod_twintau0_PerTwinFamily = 0.0_pReal
constitutive_titanmod_twingamma0_PerTwinFamily = 0.0_pReal
constitutive_titanmod_twinLambdaSlipPerTwinFamily = 0.0_pReal
constitutive_titanmod_twinp_PerTwinFamily = 0.0_pReal
constitutive_titanmod_twinq_PerTwinFamily = 0.0_pReal
constitutive_titanmod_twinq_PerTwinFamily = 0.0_pReal
allocate(constitutive_titanmod_interactionSlipSlip(lattice_maxNinteraction,maxNinstance))
constitutive_titanmod_interactionSlipSlip = 0.0_pReal
allocate(constitutive_titanmod_interaction_ee(lattice_maxNinteraction,maxNinstance))
constitutive_titanmod_interaction_ee = 0.0_pReal
allocate(constitutive_titanmod_interaction_ss(lattice_maxNinteraction,maxNinstance))
constitutive_titanmod_interaction_ss = 0.0_pReal
allocate(constitutive_titanmod_interaction_es(lattice_maxNinteraction,maxNinstance))
constitutive_titanmod_interaction_ss = 0.0_pReal
allocate(constitutive_titanmod_interactionSlipTwin(lattice_maxNinteraction,maxNinstance))
constitutive_titanmod_interactionSlipTwin = 0.0_pReal
allocate(constitutive_titanmod_interactionTwinSlip(lattice_maxNinteraction,maxNinstance))
constitutive_titanmod_interactionTwinSlip = 0.0_pReal
allocate(constitutive_titanmod_interactionTwinTwin(lattice_maxNinteraction,maxNinstance))
constitutive_titanmod_interactionSlipSlip = 0.0_pReal
constitutive_titanmod_interaction_ee = 0.0_pReal
constitutive_titanmod_interaction_ss = 0.0_pReal
constitutive_titanmod_interaction_ss = 0.0_pReal
constitutive_titanmod_interactionSlipTwin = 0.0_pReal
constitutive_titanmod_interactionTwinSlip = 0.0_pReal
constitutive_titanmod_interactionTwinTwin = 0.0_pReal
constitutive_titanmod_interactionTwinTwin = 0.0_pReal
!* Readout data from material.config file
rewind(file)
line = ''
section = 0
section = 0_pInt
write(6,*) 'titanmod: Reading material parameters from material config file'
@ -447,7 +446,7 @@ enddo
constitutive_titanmod_rho_edge0(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_rho_edge0(1:4,i)
case ('rho_screw0')
forall (j = 1:lattice_maxNslipFamily) &
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_rho_screw0(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_rho_screw0(1:4,i)
case ('slipburgers')
@ -606,7 +605,7 @@ write(6,*) 'Material Property reading done'
if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt,e=i)
if (sum(constitutive_titanmod_Nslip(:,i)) <= 0_pInt) call IO_error(231_pInt,e=i,ext_msg='nslip')
if (sum(constitutive_titanmod_Ntwin(:,i)) < 0_pInt) call IO_error(231_pInt,e=i,ext_msg='ntwin')
do f = 1,lattice_maxNslipFamily
do f = 1_pInt,lattice_maxNslipFamily
if (constitutive_titanmod_Nslip(f,i) > 0_pInt) then
if (constitutive_titanmod_rho_edge0(f,i) < 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='rho_edge0')
if (constitutive_titanmod_rho_screw0(f,i) < 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='rho_screw0')
@ -622,7 +621,7 @@ write(6,*) 'Material Property reading done'
call IO_error(231_pInt,e=i,ext_msg='kinkCriticalLength')
endif
enddo
do f = 1,lattice_maxNtwinFamily
do f = 1_pInt,lattice_maxNtwinFamily
if (constitutive_titanmod_Ntwin(f,i) > 0_pInt) then
if (constitutive_titanmod_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinburgers')
if (constitutive_titanmod_twinf0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinf0')
@ -765,7 +764,7 @@ do i = 1_pInt,maxNinstance
write(6,*) 'Determined size of state and dot state'
!* Determine size of postResults array
do o = 1,constitutive_titanmod_Noutput(i)
do o = 1_pInt,constitutive_titanmod_Noutput(i)
select case(constitutive_titanmod_output(o,i))
case('rhoedge', &
'rhoscrew', &
@ -849,20 +848,20 @@ write(6,*) 'Determining elasticity matrix'
do k=1_pInt,constitutive_titanmod_Ntwin(j,i)
do l=1_pInt,3_pInt ; do m=1_pInt,3_pInt ; do n=1_pInt,3_pInt ; do o=1_pInt,3_pInt
do p=1_pInt,3_pInt ; do q=1_pInt,3_pInt ; do r=1_pInt,3_pInt ; do s=1_pInt,3_pInt
constitutive_titanmod_Ctwin_3333(l,m,n,o,sum(constitutive_titanmod_Nslip(1:j-1,i))+k,i) = &
constitutive_titanmod_Ctwin_3333(l,m,n,o,sum(constitutive_titanmod_Nslip(1:j-1,i))+k,i) + &
constitutive_titanmod_Ctwin_3333(l,m,n,o,sum(constitutive_titanmod_Nslip(1:j-1_pInt,i))+k,i) = &
constitutive_titanmod_Ctwin_3333(l,m,n,o,sum(constitutive_titanmod_Nslip(1:j-1_pInt,i))+k,i) + &
constitutive_titanmod_Cslip_3333(p,q,r,s,i)*&
lattice_Qtwin(l,p,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* &
lattice_Qtwin(m,q,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* &
lattice_Qtwin(n,r,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* &
lattice_Qtwin(o,s,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)
lattice_Qtwin(l,p,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure)* &
lattice_Qtwin(m,q,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure)* &
lattice_Qtwin(n,r,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure)* &
lattice_Qtwin(o,s,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure)
enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo
constitutive_titanmod_Ctwin_66(:,:,k,i) = math_Mandel3333to66(constitutive_titanmod_Ctwin_3333(:,:,:,:,k,i))
enddo
enddo
!* Burgers vector, dislocation velocity prefactor for each slip system
do s = 1,constitutive_titanmod_totalNslip(i)
do s = 1_pInt,constitutive_titanmod_totalNslip(i)
f = constitutive_titanmod_slipFamily(s,i)
constitutive_titanmod_burgersPerSlipSystem(s,i) = constitutive_titanmod_burgersPerSlipFamily(f,i)
constitutive_titanmod_f0_PerSlipSystem(s,i) = constitutive_titanmod_f0_PerSlipFamily(f,i)
@ -882,7 +881,7 @@ write(6,*) 'Determining elasticity matrix'
enddo
!* Burgers vector, nucleation rate prefactor and twin size for each twin system
do t = 1,constitutive_titanmod_totalNtwin(i)
do t = 1_pInt,constitutive_titanmod_totalNtwin(i)
f = constitutive_titanmod_twinFamily(t,i)
constitutive_titanmod_burgersPerTwinSystem(t,i) = constitutive_titanmod_burgersPerTwinFamily(f,i)
constitutive_titanmod_twinsizePerTwinSystem(t,i) = constitutive_titanmod_twinsizePerTwinFamily(f,i)
@ -897,56 +896,44 @@ write(6,*) 'Determining elasticity matrix'
enddo
!* Construction of interaction matrices
do s1 = 1,constitutive_titanmod_totalNslip(i)
do s2 = 1,constitutive_titanmod_totalNslip(i)
do s1 = 1_pInt,constitutive_titanmod_totalNslip(i)
do s2 = 1_pInt,constitutive_titanmod_totalNslip(i)
constitutive_titanmod_interactionMatrixSlipSlip(s1,s2,i) = &
constitutive_titanmod_interactionSlipSlip(lattice_interactionSlipSlip(constitutive_titanmod_slipSystemLattice(s1,i), &
constitutive_titanmod_slipSystemLattice(s2,i), &
myStructure),i)
enddo; enddo
do s1 = 1,constitutive_titanmod_totalNslip(i)
do s2 = 1,constitutive_titanmod_totalNslip(i)
constitutive_titanmod_interactionMatrix_ee(s1,s2,i) = &
constitutive_titanmod_interaction_ee(lattice_interactionSlipSlip(constitutive_titanmod_slipSystemLattice(s1,i), &
constitutive_titanmod_slipSystemLattice(s2,i), &
myStructure),i)
enddo; enddo
do s1 = 1,constitutive_titanmod_totalNslip(i)
do s2 = 1,constitutive_titanmod_totalNslip(i)
constitutive_titanmod_interactionMatrix_ss(s1,s2,i) = &
constitutive_titanmod_interaction_ss(lattice_interactionSlipSlip(constitutive_titanmod_slipSystemLattice(s1,i), &
constitutive_titanmod_slipSystemLattice(s2,i), &
myStructure),i)
enddo; enddo
do s1 = 1,constitutive_titanmod_totalNslip(i)
do s2 = 1,constitutive_titanmod_totalNslip(i)
constitutive_titanmod_interactionMatrix_es(s1,s2,i) = &
constitutive_titanmod_interaction_es(lattice_interactionSlipSlip(constitutive_titanmod_slipSystemLattice(s1,i), &
constitutive_titanmod_slipSystemLattice(s2,i), &
myStructure),i)
enddo; enddo
do s1 = 1,constitutive_titanmod_totalNslip(i)
do t2 = 1,constitutive_titanmod_totalNtwin(i)
do s1 = 1_pInt,constitutive_titanmod_totalNslip(i)
do t2 = 1_pInt,constitutive_titanmod_totalNtwin(i)
constitutive_titanmod_interactionMatrixSlipTwin(s1,t2,i) = &
constitutive_titanmod_interactionSlipTwin(lattice_interactionSlipTwin(constitutive_titanmod_slipSystemLattice(s1,i), &
constitutive_titanmod_twinSystemLattice(t2,i), &
myStructure),i)
enddo; enddo
do t1 = 1,constitutive_titanmod_totalNtwin(i)
do s2 = 1,constitutive_titanmod_totalNslip(i)
do t1 = 1_pInt,constitutive_titanmod_totalNtwin(i)
do s2 = 1_pInt,constitutive_titanmod_totalNslip(i)
constitutive_titanmod_interactionMatrixTwinSlip(t1,s2,i) = &
constitutive_titanmod_interactionTwinSlip(lattice_interactionTwinSlip(constitutive_titanmod_twinSystemLattice(t1,i), &
constitutive_titanmod_slipSystemLattice(s2,i), &
myStructure),i)
enddo; enddo
do t1 = 1,constitutive_titanmod_totalNtwin(i)
do t2 = 1,constitutive_titanmod_totalNtwin(i)
do t1 = 1_pInt,constitutive_titanmod_totalNtwin(i)
do t2 = 1_pInt,constitutive_titanmod_totalNtwin(i)
constitutive_titanmod_interactionMatrixTwinTwin(t1,t2,i) = &
constitutive_titanmod_interactionTwinTwin(lattice_interactionTwinTwin(constitutive_titanmod_twinSystemLattice(t1,i), &
constitutive_titanmod_twinSystemLattice(t2,i), &
@ -954,16 +941,12 @@ write(6,*) 'Determining elasticity matrix'
enddo; enddo
!* Calculation of forest projections for edge dislocations
do s1 = 1,constitutive_titanmod_totalNslip(i)
do s2 = 1,constitutive_titanmod_totalNslip(i)
do s1 = 1_pInt,constitutive_titanmod_totalNslip(i)
do s2 = 1_pInt,constitutive_titanmod_totalNslip(i)
constitutive_titanmod_forestProjectionEdge(s1,s2,i) = &
abs(math_mul3x3(lattice_sn(:,constitutive_titanmod_slipSystemLattice(s1,i),myStructure), &
lattice_st(:,constitutive_titanmod_slipSystemLattice(s2,i),myStructure)))
enddo; enddo
!* Calculation of forest projections for screw dislocations
do s1 = 1,constitutive_titanmod_totalNslip(i)
do s2 = 1,constitutive_titanmod_totalNslip(i)
constitutive_titanmod_forestProjectionScrew(s1,s2,i) = &
abs(math_mul3x3(lattice_sn(:,constitutive_titanmod_slipSystemLattice(s1,i),myStructure), &
lattice_sd(:,constitutive_titanmod_slipSystemLattice(s2,i),myStructure)))
@ -971,22 +954,17 @@ write(6,*) 'Determining elasticity matrix'
!* Calculation of forest projections for edge dislocations in twin system
do t1 = 1,constitutive_titanmod_totalNtwin(i)
do t2 = 1,constitutive_titanmod_totalNtwin(i)
do t1 = 1_pInt,constitutive_titanmod_totalNtwin(i)
do t2 = 1_pInt,constitutive_titanmod_totalNtwin(i)
constitutive_titanmod_TwinforestProjectionEdge(t1,t2,i) = &
abs(math_mul3x3(lattice_tn(:,constitutive_titanmod_twinSystemLattice(t1,i),myStructure), &
lattice_tt(:,constitutive_titanmod_twinSystemLattice(t2,i),myStructure)))
enddo; enddo
!* Calculation of forest projections for screw dislocations in twin system
do t1 = 1,constitutive_titanmod_totalNtwin(i)
do t2 = 1,constitutive_titanmod_totalNtwin(i)
constitutive_titanmod_TwinforestProjectionScrew(t1,t2,i) = &
abs(math_mul3x3(lattice_tn(:,constitutive_titanmod_twinSystemLattice(t1,i),myStructure), &
lattice_td(:,constitutive_titanmod_twinSystemLattice(t2,i),myStructure)))
enddo; enddo
enddo
write(6,*) 'Init All done'
return
@ -998,7 +976,6 @@ function constitutive_titanmod_stateInit(myInstance)
!* initial microstructural state *
!*********************************************************************
use prec, only: pReal,pInt
use math, only: pi
use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily
implicit none
@ -1024,7 +1001,7 @@ constitutive_titanmod_stateInit = 0.0_pReal
!* Initialize basic slip state variables
! For slip
s1 = 0_pInt
do f = 1,lattice_maxNslipFamily
do f = 1_pInt,lattice_maxNslipFamily
s0 = s1 + 1_pInt
s1 = s0 + constitutive_titanmod_Nslip(f,myInstance) - 1_pInt
do s = s0,s1
@ -1037,7 +1014,7 @@ enddo
!* Initialize basic slip state variables
! For twin
ts1 = 0_pInt
do tf = 1,lattice_maxNtwinFamily
do tf = 1_pInt,lattice_maxNtwinFamily
ts0 = ts1 + 1_pInt
ts1 = ts0 + constitutive_titanmod_Ntwin(tf,myInstance) - 1_pInt
do ts = ts0,ts1
@ -1046,44 +1023,44 @@ do tf = 1,lattice_maxNtwinFamily
enddo
constitutive_titanmod_stateInit(1:ns) = rho_edge0
constitutive_titanmod_stateInit(ns+1:2*ns) = rho_screw0
constitutive_titanmod_stateInit(2*ns+1:3*ns) = shear_system0
constitutive_titanmod_stateInit(3*ns+1:3*ns+nt) = twingamma_dot0
constitutive_titanmod_stateInit(ns+1_pInt:2_pInt*ns) = rho_screw0
constitutive_titanmod_stateInit(2_pInt*ns+1_pInt:3_pInt*ns) = shear_system0
constitutive_titanmod_stateInit(3_pInt*ns+1_pInt:3_pInt*ns+nt) = twingamma_dot0
!* Initialize dependent slip microstructural variables
forall (s = 1:ns) &
forall (s = 1_pInt:ns) &
segment_edge0(s) = constitutive_titanmod_CeLambdaSlipPerSlipSystem(s,myInstance)/ &
sqrt(dot_product((rho_edge0),constitutive_titanmod_forestProjectionEdge(1:ns,s,myInstance))+ &
dot_product((rho_screw0),constitutive_titanmod_forestProjectionScrew(1:ns,s,myInstance)))
constitutive_titanmod_stateInit(3*ns+nt+1:4*ns+nt) = segment_edge0
constitutive_titanmod_stateInit(3_pInt*ns+nt+1_pInt:4_pInt*ns+nt) = segment_edge0
forall (s = 1:ns) &
forall (s = 1_pInt:ns) &
segment_screw0(s) = constitutive_titanmod_CsLambdaSlipPerSlipSystem(s,myInstance)/ &
sqrt(dot_product((rho_edge0),constitutive_titanmod_forestProjectionEdge(1:ns,s,myInstance))+ &
dot_product((rho_screw0),constitutive_titanmod_forestProjectionScrew(1:ns,s,myInstance)))
constitutive_titanmod_stateInit(4*ns+nt+1:5*ns+nt) = segment_screw0
constitutive_titanmod_stateInit(4_pInt*ns+nt+1_pInt:5_pInt*ns+nt) = segment_screw0
forall (s = 1:ns) &
forall (s = 1_pInt:ns) &
resistance_edge0(s) = &
constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerSlipSystem(s,myInstance)* &
sqrt(dot_product((rho_edge0),constitutive_titanmod_interactionMatrix_ee(1:ns,s,myInstance))+dot_product((rho_screw0), &
constitutive_titanmod_interactionMatrix_es(1:ns,s,myInstance)))
constitutive_titanmod_stateInit(5*ns+nt+1:6*ns+nt) = resistance_edge0
constitutive_titanmod_stateInit(5_pInt*ns+nt+1_pInt:6_pInt*ns+nt) = resistance_edge0
forall (s = 1:ns) &
forall (s = 1_pInt:ns) &
resistance_screw0(s) = &
constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerSlipSystem(s,myInstance)* &
sqrt(dot_product((rho_edge0),constitutive_titanmod_interactionMatrix_es(1:ns,s,myInstance))+dot_product((rho_screw0), &
constitutive_titanmod_interactionMatrix_ss(1:ns,s,myInstance)))
constitutive_titanmod_stateInit(6*ns+nt+1:7*ns+nt) = resistance_screw0
constitutive_titanmod_stateInit(6_pInt*ns+nt+1_pInt:7_pInt*ns+nt) = resistance_screw0
forall (t = 1:nt) &
forall (t = 1_pInt:nt) &
resistance_twin0(t) = 0.0_pReal
constitutive_titanmod_stateInit(7*ns+nt+1:7*ns+2*nt)=resistance_twin0
constitutive_titanmod_stateInit(7_pInt*ns+nt+1_pInt:7_pInt*ns+2_pInt*nt)=resistance_twin0
return
end function
@ -1133,8 +1110,8 @@ ns = constitutive_titanmod_totalNslip(myInstance)
nt = constitutive_titanmod_totalNtwin(myInstance)
!* Total twin volume fraction
do i=1,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ &
do i=1_pInt,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance)
enddo
!sumf = sum(state(g,ip,el)%p((6*ns+7*nt+1):(6*ns+8*nt))) ! safe for nt == 0
@ -1142,7 +1119,7 @@ sumf = sum(abs(volumefraction_pertwinsystem(1:nt))) ! safe for nt == 0
!* Homogenized elasticity matrix
constitutive_titanmod_homogenizedC = (1.0_pReal-sumf)*constitutive_titanmod_Cslip_66(:,:,myInstance)
do i=1,nt
do i=1_pInt,nt
constitutive_titanmod_homogenizedC = &
! constitutive_titanmod_homogenizedC + state(g,ip,el)%p(6*ns+7*nt+i)*constitutive_titanmod_Ctwin_66(:,:,i,myInstance)
constitutive_titanmod_homogenizedC + volumefraction_pertwinsystem(i)*constitutive_titanmod_Ctwin_66(:,:,i,myInstance)
@ -1163,10 +1140,8 @@ subroutine constitutive_titanmod_microstructure(Temperature,state,g,ip,el)
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use math, only: pi
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
use lattice, only: lattice_interactionSlipTwin,lattice_interactionTwinTwin
!use debug, only: debugger
implicit none
@ -1205,8 +1180,8 @@ nt = constitutive_titanmod_totalNtwin(myInstance)
!* State: 13*ns+2*nt+1 : 14*ns+2*nt StressRatio_screw_p
!* Total twin volume fraction
do i=1,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ &
do i=1_pInt,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance)
enddo
@ -1224,45 +1199,45 @@ sfe = 0.0002_pReal*Temperature-0.0396_pReal
! state(g,ip,el)%p(2*ns+t)/constitutive_titanmod_twinsizePerTwinSystem(t,myInstance)
! average segment length for edge dislocations in matrix
forall (s = 1:ns) &
state(g,ip,el)%p(3*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSystem(s,myInstance)/ &
forall (s = 1_pInt:ns) &
state(g,ip,el)%p(3_pInt*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSystem(s,myInstance)/ &
sqrt(dot_product(state(g,ip,el)%p(1:ns), &
constitutive_titanmod_forestProjectionEdge(1:ns,s,myInstance))+ &
dot_product(state(g,ip,el)%p(ns+1:2*ns), &
dot_product(state(g,ip,el)%p(ns+1_pInt:2_pInt*ns), &
constitutive_titanmod_forestProjectionScrew(1:ns,s,myInstance)))
! average segment length for screw dislocations in matrix
forall (s = 1:ns) &
state(g,ip,el)%p(4*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSystem(s,myInstance)/ &
forall (s = 1_pInt:ns) &
state(g,ip,el)%p(4_pInt*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSystem(s,myInstance)/ &
sqrt(dot_product(state(g,ip,el)%p(1:ns), &
constitutive_titanmod_forestProjectionEdge(1:ns,s,myInstance))+ &
dot_product(state(g,ip,el)%p(ns+1:2*ns), &
dot_product(state(g,ip,el)%p(ns+1_pInt:2_pInt*ns), &
constitutive_titanmod_forestProjectionScrew(1:ns,s,myInstance)))
!* threshold stress or slip resistance for edge dislocation motion
forall (s = 1:ns) &
state(g,ip,el)%p(5*ns+nt+s) = &
forall (s = 1_pInt:ns) &
state(g,ip,el)%p(5_pInt*ns+nt+s) = &
constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerSlipSystem(s,myInstance)*&
sqrt(dot_product((state(g,ip,el)%p(1:ns)),&
constitutive_titanmod_interactionMatrix_ee(1:ns,s,myInstance))+ &
dot_product((state(g,ip,el)%p(ns+1:2*ns)),&
dot_product((state(g,ip,el)%p(ns+1_pInt:2_pInt*ns)),&
constitutive_titanmod_interactionMatrix_es(1:ns,s,myInstance)))
!* threshold stress or slip resistance for screw dislocation motion
forall (s = 1:ns) &
state(g,ip,el)%p(6*ns+nt+s) = &
forall (s = 1_pInt:ns) &
state(g,ip,el)%p(6_pInt*ns+nt+s) = &
constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerSlipSystem(s,myInstance)*&
sqrt(dot_product((state(g,ip,el)%p(1:ns)),&
constitutive_titanmod_interactionMatrix_es(1:ns,s,myInstance))+ &
dot_product((state(g,ip,el)%p(ns+1:2*ns)),&
dot_product((state(g,ip,el)%p(ns+1_pInt:2_pInt*ns)),&
constitutive_titanmod_interactionMatrix_ss(1:ns,s,myInstance)))
!* threshold stress or slip resistance for dislocation motion in twin
forall (t = 1:nt) &
state(g,ip,el)%p(7*ns+nt+t) = &
forall (t = 1_pInt:nt) &
state(g,ip,el)%p(7_pInt*ns+nt+t) = &
constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerTwinSystem(t,myInstance)*&
(dot_product((abs(state(g,ip,el)%p(2*ns+1:2*ns+nt))),&
(dot_product((abs(state(g,ip,el)%p(2_pInt*ns+1_pInt:2_pInt*ns+nt))),&
constitutive_titanmod_interactionMatrixTwinTwin(1:nt,t,myInstance)))
@ -1288,8 +1263,8 @@ use prec, only: pReal,pInt,p_vec
use math, only: math_Plain3333to99
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem, lattice_Stwin
implicit none
!* Input-Output variables
@ -1317,8 +1292,8 @@ myStructure = constitutive_titanmod_structure(myInstance)
ns = constitutive_titanmod_totalNslip(myInstance)
nt = constitutive_titanmod_totalNtwin(myInstance)
do i=1,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ &
do i=1_pInt,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance)
enddo
@ -1336,9 +1311,9 @@ gdot_slip_edge = 0.0_pReal
gdot_slip_screw = 0.0_pReal
dgdot_dtauslip = 0.0_pReal
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_titanmod_Nslip(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_titanmod_Nslip(f,myInstance) ! process each (active) slip system in family
j = j+1_pInt
!* Calculation of Lp
@ -1347,9 +1322,9 @@ do f = 1,lattice_maxNslipFamily ! loop over all
!*************************************************
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if(myStructure>=3.and.j>3) then ! for all non-basal slip systems
if(myStructure==3) then ! only for prismatic and pyr <a> systems in hex
if(myStructure==3_pInt) then ! only for prismatic and pyr <a> systems in hex
screwvelocity_prefactor=constitutive_titanmod_debyefrequency(myInstance)* &
state(g,ip,el)%p(4*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSystem(j,myInstance)/ &
state(g,ip,el)%p(4_pInt*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSystem(j,myInstance)/ &
constitutive_titanmod_kinkcriticallength_PerSlipSystem(j,myInstance))**2
!* Stress ratio for screw ! No slip resistance for screw dislocations, only Peierls stress
@ -1472,7 +1447,7 @@ do f = 1,lattice_maxNslipFamily ! loop over all
Lp = Lp + (1.0_pReal - sumf)*gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,myStructure)
!* Calculation of the tangent of Lp
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = &
dLp_dTstar3333(k,l,m,n) + dgdot_dtauslip(j)*&
lattice_Sslip(k,l,index_myFamily+i,myStructure)*&
@ -1484,9 +1459,9 @@ enddo
gdot_twin = 0.0_pReal
dgdot_dtautwin = 0.0_pReal
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all slip families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_titanmod_Ntwin(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all slip families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_titanmod_Ntwin(f,myInstance) ! process each (active) slip system in family
j = j+1_pInt
!* Calculation of Lp
@ -1553,7 +1528,7 @@ do f = 1,lattice_maxNtwinFamily ! loop over all
Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,myStructure)
!* Calculation of the tangent of Lp
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = &
dLp_dTstar3333(k,l,m,n) + dgdot_dtautwin(j)*&
lattice_Stwin(k,l,index_myFamily+i,myStructure)*&
@ -1592,11 +1567,10 @@ function constitutive_titanmod_dotState(Tstar_v,Temperature,state,g,ip,el)
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use math, only: pi
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin
use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem, lattice_Stwin_v
implicit none
!* Input-Output variables
@ -1623,20 +1597,20 @@ MyStructure = constitutive_titanmod_structure(myInstance)
ns = constitutive_titanmod_totalNslip(myInstance)
nt = constitutive_titanmod_totalNtwin(myInstance)
do i=1,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ &
do i=1_pInt,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance)
enddo
sumf = sum(abs(volumefraction_pertwinsystem(1:nt))) ! safe for nt == 0
sumf = sum(abs(volumefraction_pertwinsystem(1_pInt:nt))) ! safe for nt == 0
constitutive_titanmod_dotState = 0.0_pReal
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_titanmod_Nslip(f,myInstance) ! process each (active) slip system in family
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_titanmod_Nslip(f,myInstance) ! process each (active) slip system in family
j = j+1_pInt
!* Multiplication of edge dislocations
@ -1667,9 +1641,9 @@ constitutive_titanmod_dotState = 0.0_pReal
!* Twin fraction evolution
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,MyStructure)) ! at which index starts my family
do i = 1,constitutive_titanmod_Ntwin(f,myInstance) ! process each (active) twin system in family
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,MyStructure)) ! at which index starts my family
do i = 1_pInt,constitutive_titanmod_Ntwin(f,myInstance) ! process each (active) twin system in family
j = j+1_pInt
!*************************************************************************
@ -1769,11 +1743,8 @@ pure function constitutive_titanmod_postResults(Tstar_v,Temperature,dt,state,g,i
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use math, only: pi
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput
use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin
implicit none
!* Definition of variables
@ -1794,8 +1765,8 @@ myStructure = constitutive_titanmod_structure(myInstance)
ns = constitutive_titanmod_totalNslip(myInstance)
nt = constitutive_titanmod_totalNtwin(myInstance)
do i=1,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ &
do i=1_pInt,nt
volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ &
constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance)
enddo
@ -1809,99 +1780,99 @@ sumf = sum(abs(volumefraction_pertwinsystem(1:nt))) ! safe for nt == 0
c = 0_pInt
constitutive_titanmod_postResults = 0.0_pReal
do o = 1,phase_Noutput(material_phase(g,ip,el))
do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
select case(constitutive_titanmod_output(o,myInstance))
case ('rhoedge')
constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p(1_pInt:ns)
c = c + ns
case ('rhoscrew')
constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p(ns+1:2*ns)
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p(ns+1_pInt:2_pInt*ns)
c = c + ns
case ('segment_edge')
constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((3_pInt*ns+nt+1_pInt):(4_pInt*ns+nt))
c = c + ns
case ('segment_screw')
constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((4*ns+nt+1):(5*ns+nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((4_pInt*ns+nt+1_pInt):(5_pInt*ns+nt))
c = c + ns
case ('resistance_edge')
constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((5*ns+nt+1):(6*ns+nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((5_pInt*ns+nt+1_pInt):(6_pInt*ns+nt))
c = c + ns
case ('resistance_screw')
constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((6*ns+nt+1):(7*ns+nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((6_pInt*ns+nt+1_pInt):(7_pInt*ns+nt))
c = c + ns
case ('velocity_edge')
constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((7*ns+2*nt+1):(8*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((7*ns+2*nt+1):(8*ns+2*nt))
c = c + ns
case ('velocity_screw')
constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((8*ns+2*nt+1):(9*ns+2*nt))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((8*ns+2*nt+1):(9*ns+2*nt))
c = c + ns
case ('tau_slip')
constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((9*ns+2*nt+1):(10*ns+2*nt)))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((9*ns+2*nt+1):(10*ns+2*nt)))
c = c + ns
case ('gdot_slip_edge')
constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((10*ns+2*nt+1):(11*ns+2*nt)))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((10*ns+2*nt+1):(11*ns+2*nt)))
c = c + ns
case ('gdot_slip_screw')
constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((11*ns+2*nt+1):(12*ns+2*nt)))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((11*ns+2*nt+1):(12*ns+2*nt)))
c = c + ns
case ('gdot_slip')
constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((10*ns+2*nt+1):(11*ns+2*nt))) + &
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((10*ns+2*nt+1):(11*ns+2*nt))) + &
abs(state(g,ip,el)%p((11*ns+2*nt+1):(12*ns+2*nt)))
c = c + ns
case ('stressratio_edge_p')
constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((12*ns+2*nt+1):(13*ns+2*nt)))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((12*ns+2*nt+1):(13*ns+2*nt)))
c = c + ns
case ('stressratio_screw_p')
constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((13*ns+2*nt+1):(14*ns+2*nt)))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((13*ns+2*nt+1):(14*ns+2*nt)))
c = c + ns
case ('shear_system')
constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((2*ns+1):(3*ns)))
constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((2*ns+1):(3*ns)))
c = c + ns
case ('shear_basal')
constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+1):(2*ns+3))))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+1):(2*ns+3))))
c = c + 1_pInt
case ('shear_prism')
constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+4):(2*ns+6))))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+4):(2*ns+6))))
c = c + 1_pInt
case ('shear_pyra')
constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+7):(2*ns+12))))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+7):(2*ns+12))))
c = c + 1_pInt
case ('shear_pyrca')
constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+13):(2*ns+24))))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+13):(2*ns+24))))
c = c + 1_pInt
case ('rhoedge_basal')
constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((1):(3)))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((1):(3)))
c = c + 1_pInt
case ('rhoedge_prism')
constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((4):(6)))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((4):(6)))
c = c + 1_pInt
case ('rhoedge_pyra')
constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((7):(12)))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((7):(12)))
c = c + 1_pInt
case ('rhoedge_pyrca')
constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((13):(24)))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((13):(24)))
c = c + 1_pInt
case ('rhoscrew_basal')
constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((ns+1):(ns+3)))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((ns+1):(ns+3)))
c = c + 1_pInt
case ('rhoscrew_prism')
constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((ns+4):(ns+6)))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((ns+4):(ns+6)))
c = c + 1_pInt
case ('rhoscrew_pyra')
constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((ns+7):(ns+12)))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((ns+7):(ns+12)))
c = c + 1_pInt
case ('rhoscrew_pyrca')
constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((ns+13):(ns+24)))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((ns+13):(ns+24)))
c = c + 1_pInt
case ('shear_total')
constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+1):(3*ns))))
c = c + 1
constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+1):(3*ns))))
c = c + 1_pInt
case ('twin_fraction')
constitutive_titanmod_postResults(c+1:c+nt) = abs(volumefraction_pertwinsystem(1:nt))
constitutive_titanmod_postResults(c+1_pInt:c+nt) = abs(volumefraction_pertwinsystem(1:nt))
c = c + nt
! 'rhoedge', &
! 'rhoscrew', &