nucleation at boundary

This commit is contained in:
achalhp 2024-06-10 14:03:35 +05:30
parent 872228dd47
commit b26e60b28f
1 changed files with 88 additions and 29 deletions

View File

@ -287,6 +287,7 @@ do ph = 1, phases%length
allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers)) allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers))
allocate(geom(ph)%IParea(nIPneighbors,Nmembers)) allocate(geom(ph)%IParea(nIPneighbors,Nmembers))
allocate(geom(ph)%IPcoordinates(3,Nmembers)) allocate(geom(ph)%IPcoordinates(3,Nmembers))
call storeGeometry(ph)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
@ -500,16 +501,16 @@ real(pReal), dimension(3,3), intent(out) :: &
integer :: & integer :: &
n, & ! neighbor index n, & ! neighbor index
neighbor_e, & ! element index of my neighbor neighbor_e, neighbor_e1, & ! element index of my neighbor
neighbor_i, & ! integration point index of my neighbor neighbor_i, & ! integration point index of my neighbor
neighbor_me, & neighbor_me, &
neighbor_phase neighbor_phase
real(pReal) :: & real(pReal) :: &
random, & random, random1, &
nRealNeighbors !Achal number of really existing neighbors nRealNeighbors !Achal number of really existing neighbors
integer :: & integer :: &
twin_var twin_var, var_growth
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
fdot_twin fdot_twin
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
@ -520,19 +521,10 @@ deltaFp = math_I3
!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities !Achal !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities !Achal
associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dlt => deltastate(ph)) associate(prm => param(ph), stt => state(ph), dlt => deltastate(ph))
nRealNeighbors = 0.0_pReal !Achal nRealNeighbors = 0.0_pReal !Achal
neighbors: do n = 1,nIPneighbors
neighbor_e = geom(ph)%IPneighborhood(1,n,en)
!write(6,*) 'neighbor_e', neighbor_e
neighbor_i = geom(ph)%IPneighborhood(2,n,en)
neighbor_me = material_phaseEntry(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) !Neighbour offset
neighbor_phase = material_phaseID(1,(neighbor_e-1)*discretization_nIPs + neighbor_i)
enddo neighbors
!tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)] !tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)]
!twin_var = maxloc((0.05_pReal*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw)/prm%gamma_char,dim=1) ! This prints values from 1 to 6, fdot0_twin is taken as 0.05 !twin_var = maxloc((0.05_pReal*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw)/prm%gamma_char,dim=1) ! This prints values from 1 to 6, fdot0_twin is taken as 0.05
twin_var = maxloc(stt%f_twin(:,en),dim=1) twin_var = maxloc(stt%f_twin(:,en),dim=1)
@ -544,23 +536,55 @@ associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dlt => deltas
call RANDOM_NUMBER(random) call RANDOM_NUMBER(random)
call RANDOM_NUMBER(random1)
!write(6,*)'random',random !delete this !write(6,*)'random',random !delete this
!if (en==1) write(6,*)'f_twin', stt%f_twin(:,en) !if (en==1) write(6,*)'f_twin', stt%f_twin(:,en)
Ability_Nucleation: if(stt%f_twin(twin_var,en)>(stt%fmc_twin(twin_var,en)+prm%checkstep)) then
stt%fmc_twin(twin_var,en) = stt%fmc_twin(twin_var,en)+prm%checkstep
Success_Nucleation: if (random <= stt%f_twin(twin_var,en)) then ! Instead of sum take max
twinJump = .true.
deltaFp = prm%CorrespondanceMatrix(:,:,twin_var)
write(6,*)'deltaFp',deltaFp,'element',en
dlt%f_twin(:,en) = 0.0_pReal - stt%f_twin(:,en)
dlt%fmc_twin(:,en) = 0.0_pReal - stt%fmc_twin(:,en)
dlt%frozen(en) = 1.0_pReal - stt%frozen(en)
dlt%variant_twin(en) = twin_var - stt%variant_twin(en) ! Achal LHS is real, RHS integer ! why this equation?
end if Success_Nucleation
endif Ability_Nucleation neighbors: do n = 1,nIPneighbors
neighbor_e = geom(ph)%IPneighborhood(1,n,en)
if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), &
phase_O_0(ph)%data(neighbor_e)%asQuaternion()))) then
Ability_Nucleation: if(stt%f_twin(twin_var,en)>(stt%fmc_twin(twin_var,en)+prm%checkstep)) then
stt%fmc_twin(twin_var,en) = stt%fmc_twin(twin_var,en)+prm%checkstep
Success_Nucleation: if (random <= stt%f_twin(twin_var,en)) then ! Instead of sum take max
twinJump = .true.
deltaFp = prm%CorrespondanceMatrix(:,:,twin_var)
dlt%f_twin(:,en) = 0.0_pReal - stt%f_twin(:,en)
dlt%fmc_twin(:,en) = 0.0_pReal - stt%fmc_twin(:,en)
dlt%frozen(en) = 1.0_pReal - stt%frozen(en)
dlt%variant_twin(en) = twin_var !- stt%variant_twin(en) ! Achal LHS is real, RHS integer ! why this equation?
!write(6,*)'variant_twin',stt%variant_twin(en),'element',en,'frozen',stt%frozen(en)
end if Success_Nucleation
endif Ability_Nucleation
endif
enddo neighbors
NeighborLoop: do n = 1,nIPneighbors
neighbor_e1 = geom(ph)%IPneighborhood(1,n,en)
if(stt%variant_twin(neighbor_e1)>0) then !< Check if neighbor is twinned
var_growth = stt%variant_twin(neighbor_e1)
exit NeighborLoop
endif
enddo NeighborLoop
Growth_Criteria: if(var_growth>0) then !< If neighbor twinned,
!Ability_Growth: if(stt%f_twin(twin_var,en)>(stt%fmc_twin(twin_var,en)+prm%checkstep(twin_var))) then !< Frequency control
!stt%fmc_twin(twin_var,en) = stt%fmc_twin(twin_var,en)+prm%checkstep(twin_var)
Success_Growth: if (random1 <= stt%f_twin(twin_var,en)) then !< Random sampling
twinJump = .true. !< Output flag
deltaFp = prm%CorrespondanceMatrix(:,:,twin_var) !< Correspondence Matrix
endif Success_Growth
!endif Ability_Growth
endif Growth_Criteria
end associate end associate
@ -577,12 +601,14 @@ integer, intent(in)::&
ph, & ph, &
en en
! These are updated at every strain increment. What should these initilizations be?
associate(dlt => deltastate(ph)) associate(dlt => deltastate(ph))
dlt%f_twin(:,en) = 0.0_pReal dlt%f_twin(:,en) = 0.0_pReal
dlt%fmc_twin(:,en) = 0.0_pReal dlt%fmc_twin(:,en) = 0.0_pReal
!dlt%variant_twin(en) = 1.0_pReal !dlt%variant_twin(en) = 0.0_pReal
!dlt%frozen(en) = 1.0_pReal !dlt%frozen(en) = 0.0_pReal
end associate end associate
@ -746,6 +772,7 @@ associate(prm => param(ph), stt => state(ph))
fdot_twin = (0.05_pReal*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw)/prm%gamma_char !Achal 0.05 is constant fdot_twin = (0.05_pReal*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw)/prm%gamma_char !Achal 0.05 is constant
else where else where
dot_gamma_tw = 0.0_pReal dot_gamma_tw = 0.0_pReal
fdot_twin = 0.0_pReal
end where end where
if (present(ddot_gamma_dtau_tw)) then if (present(ddot_gamma_dtau_tw)) then
@ -760,4 +787,36 @@ end associate
end subroutine kinetics_tw end subroutine kinetics_tw
subroutine storeGeometry(ph)
integer, intent(in) :: ph
integer :: ce, co, nCell
real(pReal), dimension(:), allocatable :: V
integer, dimension(:,:,:), allocatable :: neighborhood
real(pReal), dimension(:,:), allocatable :: area, coords
real(pReal), dimension(:,:,:), allocatable :: areaNormal
nCell = product(shape(IPvolume))
V = reshape(IPvolume,[nCell])
neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell])
area = reshape(IParea,[nIPneighbors,nCell])
areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell])
coords = reshape(discretization_IPcoords,[3,nCell])
do ce = 1, size(material_homogenizationEntry,1)
do co = 1, homogenization_maxNconstituents
if (material_phaseID(co,ce) == ph) then
geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
geom(ph)%IPneighborhood(:,:,material_phaseEntry(co,ce)) = neighborhood(:,:,ce)
geom(ph)%IParea(:,material_phaseEntry(co,ce)) = area(:,ce)
geom(ph)%IPareaNormal(:,:,material_phaseEntry(co,ce)) = areaNormal(:,:,ce)
geom(ph)%IPcoordinates(:,material_phaseEntry(co,ce)) = coords(:,ce)
end if
end do
end do
end subroutine
end submodule phenopowerlaw end submodule phenopowerlaw