Commit e050b342 authored by Tianle Cheng's avatar Tianle Cheng
Browse files

Update README.md

parent 7103967d
# Flood fill for phase field
# This is a FORTRAN code to identify region outlined by a phase-field variable
module flood_fill
!
implicit none
!
contains
!
!-------------------------------------------------------------------------------------
! This is a 2D code, assume the system size is
! System_Size(1)xSystem_Size(2)xSystem_Size(3)
!
Subroutine flood_fill_PSD(yita, area)
implicit none
integer, parameter :: MAXPN=9999, System_Size(3)=(/1, 128, 128/)
real :: System_Volume=real(System_Size(1)*System_Size(2)*System_Size(3))
integer i,j,k, ka, kb, jr, jl, nprep, SSeedA(MAXPN,50,3), SSeedB(MAXPN,50,3), use_stack(MAXPN), Sstack(MAXPN, 2)
integer :: js,jt, jtt, landing_P(3)
logical :: safety_valve
real :: boundvalue=0.5
real :: area(MAXPN), yita(System_Size(1), System_Size(2), System_Size(3))
integer :: Is_scanning_P(MAXPN), counted(System_Size(1), System_Size(2), System_Size(3))
SSeedA=0; SSeedB=0; Sstack=1; use_stack=0
area=0; counted=0
nprep=1; landing_P=0
i=1; j=1; k=1
safety_valve=.TRUE.
do while (k.le.System_Size(3).and.safety_valve)
if(yita(i,j,k).ge.boundvalue.and.counted(i,j,k).eq.0)then
jr=j
else
jr=j+1
if(jr.eq.System_Size(2)+1)jr=1
do while ((yita(i,jr,k).lt.boundvalue.or.counted(i,jr,k).eq.1).and.jr.ne.j)
jr=jr+1
if(jr.eq.System_Size(2)+1) jr=1
enddo
endif
If(jr.eq.j.and.(yita(i,j,k).lt.boundvalue.or.counted(i,j,k).eq.1))then ! this means no uncounted particle in this line
k=k+1; j=1 ! if no particle in this line, move up one line and continue scanning
else ! else there is a uncounted precipitate now
! write(12,*)"found a seed at ", jr, k
Is_scanning_P(nprep)=1
jl=jr;
do while (yita(i,jl,k).ge.boundvalue.and.counted(i,jl,k).ne.1)
jl=jl-1
if(jl.eq.0)jl=System_Size(2)
enddo
jl=jl+1
if(jl.eq.System_Size(2)+1)jl=1
do while (yita(i,jr,k).ge.boundvalue.and.counted(i,jr,k).ne.1)
jr=jr+1
if(jr.eq.System_Size(2)+1)jr=1
enddo
jr=jr-1
if(jr.eq.0) jr=System_Size(2)
! Now jl and jr mark the left and right bound of particle line
if(jr.ge.jl)then
counted(i, jl:jr,k)=1
area(nprep)=area(nprep)+ jr-jl+1
else
counted(i, 1:jr,k)=1; counted(i, jl: System_Size(2),k)=1
area(nprep)=area(nprep)+jr; area(nprep)=area(nprep)+System_Size(2)-jl+1
endif
! write(12,*)"integrate line", k," jl:jr ", jl, jr, "area:", area(nprep)
if(use_stack(nprep).eq.0)then
landing_P(2)=jr+1; landing_P(3)=k
if(jr+1.eq.System_Size(2)+1)then
landing_P(2)=1; landing_P(3)=k+1
endif
endif
ka=k+1
if(ka.eq.System_Size(3)+1) ka=1
js=jl
do while(js.ne.jr.and.(yita(i,js,ka).lt.boundvalue.or.counted(i,js,ka).eq.1))
js=js+1
if(js.eq.System_Size(2)+1)js=1
enddo
if((js.ne.jr.and.yita(i,js,ka).ge.boundvalue).and.counted(i,js,ka).eq.0) then
SSeedA(nprep,Sstack(nprep,1),2)= js; SSeedA(nprep,Sstack(nprep,1),3)= ka
Sstack(nprep,1)=Sstack(nprep,1)+1
jt=js+1
if(jt.eq.System_Size(2)+1)jt=1
jtt=jt+1
if(jtt.eq.System_Size(2)+1)jtt=1
if(js.ne.jr.and.jt.ne.jr)then
do while (jtt.ne.jr)
if(yita(i,jt,ka).lt.boundvalue.and.yita(i,jtt,ka).ge.boundvalue.and.counted(i,jtt,ka).eq.0)then
SSeedA(nprep,Sstack(nprep,1),2)= jtt; SSeedA(nprep,Sstack(nprep,1),3)= ka
Sstack(nprep,1)=Sstack(nprep,1)+1
endif
jt=jt+1
if(jt.eq.System_Size(2)+1)jt=1
jtt=jt+1
if(jtt.eq.System_Size(2)+1)jtt=1
enddo
endif
else
if(js.eq.jr.and.js.eq.jl.and.yita(i,js,ka).ge.boundvalue.and.counted(i,js,ka).eq.0)then
SSeedA(nprep,Sstack(nprep,1),2)= js; SSeedA(nprep,Sstack(nprep,1),3)= ka
Sstack(nprep,1)=Sstack(nprep,1)+1
endif
endif
kb=k-1
if(kb.eq.0) kb=System_Size(3)
js=jl
do while(js.ne.jr.and.(yita(i,js,kb).lt.boundvalue.or.counted(i,js,kb).eq.1))
js=js+1
if(js.gt.System_Size(2))js=1
enddo
if((js.ne.jr.and.yita(i,js,kb).ge.boundvalue).and.counted(i,js,kb).eq.0) then
SSeedB(nprep,Sstack(nprep,2),2)= js; SSeedB(nprep,Sstack(nprep,2),3)= kb
Sstack(nprep,2)=Sstack(nprep,2)+1
jt=js+1
if(jt.eq.System_Size(2)+1)jt=1
jtt=jt+1
if(jtt.eq.System_Size(2)+1)jtt=1
if(js.ne.jr.and.jt.ne.jr)then
do while (jtt.ne.jr)
if(yita(i,jt,kb).lt.boundvalue.and.yita(i,jtt,kb).ge.boundvalue.and.counted(i,jtt,kb).eq.0)then
SSeedB(nprep,Sstack(nprep,2),2)= jtt; SSeedB(nprep,Sstack(nprep,2),3)= kb
Sstack(nprep,2)=Sstack(nprep,2)+1
endif
jt=jt+1
if(jt.eq.System_Size(2)+1)jt=1
jtt=jt+1
if(jtt.eq.System_Size(2)+1)jtt=1
enddo
endif
else
if(js.eq.jr.and.js.eq.jl.and.yita(i,js,kb).ge.boundvalue.and.counted(i,js,kb).eq.0)then
SSeedB(nprep,Sstack(nprep,2),2)=js; SSeedB(nprep,Sstack(nprep,2),3)= kb
Sstack(nprep,2)=Sstack(nprep,2)+1
endif
endif
if(Sstack(nprep,1).gt.1)then ! There is particle line above
use_stack(nprep)=1
j=SSeedA(nprep,Sstack(nprep,1)-1,2); k=SSeedA(nprep,Sstack(nprep,1)-1,3)
Sstack(nprep,1)=Sstack(nprep,1)-1
else
if(Sstack(nprep,2).gt.1)then ! There is particle line below
use_stack(nprep)=1
j=SSeedB(nprep,Sstack(nprep,2)-1,2); k=SSeedB(nprep,Sstack(nprep,2)-1,3)
Sstack(nprep,2)=Sstack(nprep,2)-1
else
Is_scanning_P(nprep)=0
j=landing_P(2); k=landing_P(3)
nprep=nprep+1
endif
endif
Endif
if(area(nprep).gt.System_volume/2.or.nprep.ge.MAXPN-2.or.Sstack(nprep,1).gt.40.or.Sstack(nprep,2).gt.40)then
safety_valve=.False.
write(12,*)'Safety valve is off, NOT ready for PSD', 'n & area(n) & Stack12:', nprep, &
area(nprep), Sstack(nprep,1:2)
endif
Enddo
End Subroutine flood_fill_PSD
end module flood_fill
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment