Next: GAUSS.F Up: Listing of key Previous: DC3D.F
c ************************ ac3d.f ********************************
c BACKGROUND
c This program accepts as input a 3-d digital image, converting it
c into a complex conductor network. The conjugate gradient method
c is used to solve this finite difference representation of Laplace's
c equation for complex conductivity problems.
c Periodic boundary conditions are maintained.
c In the comments below, (USER) means that this is a section of code
c that the user might have to change for his particular problem.
c Therefore the user is encouraged to search for this string.
c PROBLEM AND VARIABLE DEFINITION
c The mathematical problem that the conjugate gradient algorithm solves
c is the minimization of the quadratic form 1/2 uAu, where
c u is the vector of voltages, and A is generated from the bond
c conductances between pixels. Nodes are thought of as being in the
c center of pixels. The minimization is constrained by maintaining an
c general direction applied electric field across the sample.
c The vectors gx,gy,gz are bond conductances, u is the voltage array,
c and gb,h, and Ah are auxiliary variables, used in subroutine dembx.
c The vector pix contains the phase labels for each pixel.
c The small vector a(i) is the volume fraction
c of the ith phase, and currx, curry, currz are the total volume-averaged
c complex currents in the x,y, and z directions.
c DIMENSIONS
c The vectors gx,gy,gz,u,gb,h,Ah,list,pix are all dimensioned
c ns2 = (nx+2)*(ny+2)*(nz+2). This number is used, rather than the
c system size nx x ny x nz, because an extra layer of pixels is
c put around the system to be able to maintain periodic boundary
c conditions (see manual, Sec. 3.3). The arrays pix and list are also
c dimensioned this way. At present the program is set up for 100 phases,
c but that can easily be changed by the user, by changing the dimensions
c of sigma, a, and be. Note that be has both dimensions equal to
c each other. The parameter nphase gives the number of phases
c being considered. The parameter ntot is the total number of phases
c possible in the program, and should be equal to the dimension
c of sigma, a, and be. All arrays are passed to subroutines in
c the call statements.
c STRONGLY RECOMMENDED: READ MANUAL BEFORE USING THE PROGRAM!!
c (USER) Change these dimensions for different system sizes. All
c dimensions in the subroutines are passed, so do not need to be
c changed. The dimensions of sigma, a, and be should be equal to
c the value of ntot.
complex gx(10648),gy(10648),u(10648),gz(10648)
complex gb(10648),h(10648),Ah(10648)
complex currx,curry,currz,sigma(100,3),be(100,100,3)
real a(100)
integer*2 pix(10648)
integer*4 list(10648)
c (USER) Unit 9 is the microstructure input file, unit 7 is the
c results output file.
open(unit=9,file='microstructure.dat')
open(unit=7,file='outputfile.out')
c (USER) real image size is nx x ny x nz
nx=20
ny=20
nz=20
c auxiliary variables involving the system size
nx1=nx+1
ny1=ny+1
nz1=nz+1
nx2=nx+2
ny2=ny+2
nz2=nz+2
L22=nx2*ny2
write(7,1111) nx,ny,nz,nx*ny*nz
1111 format(' Image is ',3i6,' No. of real sites = ',i8)
c computational image size ns2 is nx2 x ny2 x nz2
ns2=nx2*ny2*nz2
c defines the value of pi for later use
pi=4.0*atan(1.0)
c (USER) set cutoff for norm squared of gradient, gtest. gtest is
c the stopping criterion, compared to gb*gb. When gb*gb is less
c than gtest=abc*ns2, then the rms value of the gradient at a pixel
c is less than sqrt(abc).
gtest=1.0e-16*ns2
c (USER) nphase is the number of phases being considered in the problem.
c The values of pix(m) will run from 1 to nphase. ntot is the total
c number of phases possible in the program, and is the dimension of
c sigma, a, and be.
nphase=2
ntot=100
c Make list of real (interior) sites, used in subroutine dembx. The 1-d
c labelling scheme goes over all ns2 sites, so a list of the real sites
c is needed.
nlist=0
do 103 i=2,nx1
do 102 j=2,ny1
do 101 k=2,nz1
m=i+(j-1)*nx2+(k-1)*L22
nlist=nlist+1
list(nlist)=m
101 continue
102 continue
103 continue
c Compute average current in each pixel.
c (USER) npoints is the number of microstructures to use.
c nfreq is the number of frequencies to be computed.
c The program is set up assuming that the effective
c conductivity is going to be solved for at several different
c frequencies on the same microstructure.
npoints=1
do 400 micro=1,npoints
c Read in a microstructure in subroutine ppixel, and set up
c pix(m) with the appropriate phase assignments.
call ppixel(pix,nx2,ny2,nz2,a,ns2,nphase,ntot)
c output phase volume fractions
do 99 i=1,nphase
write(7,299) i,a(i)
299 format(' Phase fraction of ',i3,' = ',f12.6)
99 continue
c (USER) Set components of applied field, E = (ex,ey,ez)
ex=1.0
ey=1.0
ez=1.0
write(7,*) 'Applied field components:'
write(7,*) 'ex = ',ex,' ey = ',ey,' ez = ',ez
c Initialize the voltage distribution by putting on uniform field.
c Only do this for the first frequency considered, thereafter use the
c previous frequency's voltages as a starting point.
do 30 k=1,nz2
do 30 j=1,ny2
do 30 i=1,nx2
m=(k-1)*nx2*ny2+nx2*(j-1)+i
u(m)=-ex*i-ey*j-ez*k
30 continue
c (USER) Set how many frequencies need to be computed.
nfreq=50
c Loop over desired frequencies.
do 300 nf=1,nfreq
c (USER) Define frequency to use each time. Alter this statement to get
c different frequencies. The frequencies are in Hz, according to
c the units used for sigma.
w=10.**((nf-1)*11.4/49.-3.)
c convert to angular frequency
w=w*2.*pi
write(7,*) 'No.',nf, ' angular frequency = ',w,' radians'
c (USER) input value of complex conductivity tensor for each phase
c (diagonal only). 1,2,3 = x,y,z, respectively.
sigma(1,1)=cmplx(1.0,10.*w)
sigma(1,2)=cmplx(1.0,10.*w)
sigma(1,3)=cmplx(1.0,10.*w)
sigma(2,1)=cmplx(0.5,1.*w)
sigma(2,2)=cmplx(0.5,1.*w)
sigma(2,3)=cmplx(0.5,1.*w)
c bond() sets up conductor network in gx,gy,gz 1-d arrays
call bond(pix,gx,gy,gz,nx2,ny2,nz2,ns2,sigma,be,nphase,ntot)
c Subroutine dembx accepts gx,gy,gz and solves for the voltage field
c that minimizes the dissipated energy. As a starting point for u,
c the voltage vector, each frequency uses the voltages left over from the
c previous minimization. This can often reduce total run time dramatically,
c compared to starting with a new voltage vector each time, as in do loop
c 30 above.
call dembx(nx2,ny2,nz2,ns2,gx,gy,gz,u,ic,gb,h,Ah,list,nlist,gtest)
c find final current after voltage solution is done
call current(nx2,ny2,nz2,ns2,currx,curry,currz,u,gx,gy,gz)
write(7,*) 'Average current in x direction= ',currx
write(7,*) 'Average current in y direction= ',curry
write(7,*) 'Average current in z direction= ',currz
write(7,*) ic,' number of conjugate gradient cycles needed'
call flush(7)
300 continue
400 continue
end
c Subroutine that performs the conjugate gradient solution routine to
c find the correct set of nodal voltages
subroutine dembx(nx2,ny2,nz2,ns2,gx,gy,gz,u,ic,gb,h,Ah,
& list,nlist,gtest)
complex gx(ns2),gy(ns2),u(ns2),gb(ns2)
complex Ah(ns2),h(ns2),gz(ns2)
complex gg,hAh,lambda,gglast,gamma,ravg,currx,curry,currz
integer*4 list(ns2),ncheck
c Note: voltage gradients are maintained because in the conjugate gradient
c relaxation algorithm, the voltage vector is only modified by adding a
c periodic vector to it.
L22=nx2*ny2
c First stage, compute initial value of gradient (gb), initialize h, the
c conjugate gradient direction, and compute norm squared of gradient vector.
call prod(nx2,ny2,nz2,ns2,gx,gy,gz,u,gb)
do 20 i=1,ns2
h(i)=gb(i)
20 continue
c Complex variable gg is the norm squared of the gradient vector
gg=cmplx(0.0,0.0)
do 105 k=1,nlist
m=list(k)
gg=gb(m)*gb(m)+gg
105 continue
c Second stage, initialize Ah variable, compute parameter lamdba,
c make first change in voltage array, update gradient (gb) vector
if(abs(real(gg)).lt.gtest) goto 44
call prod(nx2,ny2,nz2,ns2,gx,gy,gz,h,Ah)
hAh=cmplx(0.0,0.0)
do 205 k=1,nlist
m=list(k)
hAh=hAh+h(m)*Ah(m)
205 continue
lambda=gg/hAh
do 50 i=1,ns2
u(i)=u(i)-lambda*h(i)
gb(i)=gb(i)-lambda*Ah(i)
50 continue
c third stage: iterate conjugate gradient solution process until
c real(gg) < gtest criterion is satisfied.
c (USER) The parameter ncgsteps is the total number of conjugate gradient steps
c to go through. Only in very unusual problems, like when the conductivity
c of one phase is much higher than all the rest, will this many steps be
c used.
ncgsteps=30000
do 33 icc=1,ncgsteps
gglast=gg
gg=cmplx(0.0,0.0)
do 305 k=1,nlist
m=list(k)
gg=gb(m)*gb(m)+gg
305 continue
call flush(7)
if(abs(real(gg)).lt.gtest) goto 44
gamma=gg/gglast
c update conjugate gradient direction
do 70 i=1,ns2
h(i)=gb(i)+gamma*h(i)
70 continue
call prod(nx2,ny2,nz2,ns2,gx,gy,gz,h,Ah)
hAh=cmplx(0.0,0.0)
do 401 k=1,nlist
m=list(k)
hAh=hAh+h(m)*Ah(m)
401 continue
lambda=gg/hAh
c update voltage, gradient vectors
do 90 i=1,ns2
u(i)=u(i)-lambda*h(i)
gb(i)=gb(i)-lambda*Ah(i)
90 continue
c (USER) This piece of code forces dembx to write out the total current and
c the norm of the gradient squared, every ncheck conjugate gradient steps,
c in order to see how the relaxation is proceeding. If the currents become
c unchanging before the relaxation is done, then gtest was picked to be
c smaller than was necessary.
ncheck=30
if(ncheck*(icc/ncheck).eq.icc) then
write(7,*) icc
write(7,*) ' gg = ',gg
c call current subroutine
call current(nx2,ny2,nz2,ns2,currx,curry,currz,u,gx,gy,gz)
write(7,*) ' currx = ',currx
write(7,*) ' curry = ',curry
write(7,*) ' currz = ',currz
end if
call flush(7)
33 continue
write(7,*) ' Iteration failed to converge after',ncgsteps,' steps'
44 continue
ic=icc
return
end
c The matrix product subroutine
subroutine prod(nx2,ny2,nz2,ns2,gx,gy,gz,xw,yw)
complex gx(ns2),gy(ns2),xw(ns2),gz(ns2),yw(ns2)
c xw is the input vector, yw = (A)(xw) is the output vector
c auxiliary variables involving the system size
nx1=nx2-1
ny1=ny2-1
nz1=nz2-1
nx=nx2-2
ny=ny2-2
nz=nz2-2
L22=nx2*ny2
c Perform basic matrix multiplication, results in incorrect information at
c periodic boundaries.
do 10 i=1,ns2
yw(i)=cmplx(0.0,0.0)
10 continue
do 20 i=L22+1,ns2-L22
yw(i)=-xw(i)*(gx(i-1)+gx(i)+gz(i-L22)+gz(i)+gy(i)+gy(i-nx2))
yw(i)=yw(i)+gx(i-1)*xw(i-1)+gx(i)*xw(i+1)
+ +gz(i-L22)*xw(i-L22)+gz(i)*xw(i+L22)+gy(i)*xw(i+nx2)
+ +gy(i-nx2)*xw(i-nx2)
20 continue
c Correct terms at periodic boundaries (Section 3.3 in manual)
c x faces
do 30 k=1,nz2
do 30 j=1,ny2
yw((k-1)*L22+nx2*(j-1)+nx2)=yw((k-1)*L22+nx2*(j-1)+2)
yw((k-1)*L22+nx2*(j-1)+1)=yw((k-1)*L22+nx2*(j-1)+nx1)
30 continue
c y faces
do 40 k=1,nz2
do 40 i=1,nx2
yw((k-1)*L22+i)=yw((k-1)*L22+ny*nx2+i)
yw((k-1)*L22+ny1*nx2+i)=yw((k-1)*L22+nx2+i)
40 continue
c z faces
do 50 m=1,L22
yw(m)=yw(m+nz*L22)
yw(m+nz1*L22)=yw(m+L22)
50 continue
return
end
c Subroutine that determines the correct bond conductances that are used
c to compute multiplication by the matrix A
subroutine bond(pix,gx,gy,gz,nx2,ny2,nz2,ns2,sigma,be,nphase,ntot)
complex gx(ns2),gy(ns2),gz(ns2),sigma(ntot,3),be(ntot,ntot,3)
integer*2 pix(ns2)
c auxiliary variables involving the system size
L22=nx2*ny2
nx=nx2-2
ny=ny2-2
nz=nz2-2
nx1=nx2-1
ny1=ny2-1
nz1=nz2-1
c Set values of conductor for phase(i)--phase(j) interface,
c store in array be(i,j,3). If either phase i or j has zero conductivity,
c then be(i,j,3)=0.
do 10 m=1,3
do 10 i=1,nphase
do 10 j=1,nphase
if(real(sigma(i,m)).eq.0.0.and.aimag(sigma(i,m)).eq.0.0) then
be(i,j,m)=cmplx(0.0,0.0)
goto 10
end if
if(real(sigma(j,m)).eq.0.0.and.aimag(sigma(j,m)).eq.0.0) then
be(i,j,m)=cmplx(0.0,0.0)
goto 10
end if
be(i,j,m)=1./(0.5/sigma(i,m)+0.5/sigma(j,m))
10 continue
c Trim off x and y faces so that no current can flow past periodic
c boundaries. This step is not really necessary, as the voltages on the
c periodic boundaries will be matched to the corresponding real voltages
c in each conjugate gradient step.
do 20 k=1,nz2
do 15 j=1,ny2
gx((k-1)*L22+nx2*(j-1)+nx2)=cmplx(0.0,0.0)
15 continue
do 16 i=1,nx2
gy((k-1)*L22+ny1*nx2+i)=cmplx(0.0,0.0)
16 continue
20 continue
c Set up conductor network
c bulk--gz
do 30 i=1,nx2
do 30 j=1,ny2
do 30 k=1,nz1
m=(k-1)*L22+(j-1)*nx2+i
i1=i
j1=j
k1=k+1
m1=(k1-1)*L22+(j1-1)*nx2+i1
gz(m)=be(pix(m),pix(m1),3)
30 continue
c bulk---gy
do 40 i=1,nx2
do 40 j=1,ny1
do 40 k=2,nz1
m=(k-1)*L22+(j-1)*nx2+i
j1=j+1
i1=i
k1=k
m1=(k1-1)*L22+(j1-1)*nx2+i1
gy(m)=be(pix(m),pix(m1),2)
40 continue
c bulk--gx
do 50 i=1,nx1
do 50 j=1,ny2
do 50 k=2,nz1
m=(k-1)*L22+(j-1)*nx2+i
i1=i+1
j1=j
k1=k
m1=(k1-1)*L22+(j1-1)*nx2+i1
gx(m)=be(pix(m),pix(m1),1)
50 continue
return
end
c Subroutine that sets up the image, either by reading it from file,
c or generating it internally
subroutine ppixel(pix,nx2,ny2,nz2,a,ns2,nphase,ntot)
real a(ntot)
integer*2 pix(ns2)
c (USER) If you want to set up a test image inside the program, instead
c of reading it in from a file, this should be done inside this subroutine.
c auxiliary variables involving the system size
nx=nx2-2
ny=ny2-2
nz=nz2-2
L22=nx2*ny2
c Initialize phase fraction array.
do 120 i=1,nphase
a(i)=0.0
120 continue
c Use 1-d labelling scheme as shown in manual
do 100 k=2,nz2-1
do 100 j=2,ny2-1
do 100 i=2,nx2-1
m=(k-1)*L22+(j-1)*nx2+i
read(9,*) pix(m)
a(pix(m))=a(pix(m))+1.0
100 continue
do 220 i=1,nphase
a(i)=a(i)/float(nx*ny*nz)
220 continue
c now map periodic boundaries of pix (see Section 3.3, Figure 3 in manual)
do 110 k=1,nz2
do 110 j=1,ny2
do 110 i=1,nx2
if(i.gt.1.and.i.lt.nx2) then
if(j.gt.1.and.j.lt.ny2) then
if(k.gt.1.and.k.lt.nz2) then
goto 110
end if
end if
end if
k1=k
if(k.eq.1) k1=k+nz
if(k.eq.nz2) k1=k-nz
j1=j
if(j.eq.1) j1=j+ny
if(j.eq.ny2) j1=j-ny
i1=i
if(i.eq.1) i1=i+nx
if(i.eq.nx2) i1=i-nx
m=(k-1)*L22+(j-1)*nx2+i
m1=(k1-1)*L22+(j1-1)*nx2+i1
pix(m)=pix(m1)
110 continue
c Check for wrong phase labels--less than 1 or greater than nphase
do 500 m=1,ns2
if(pix(m).lt.1) then
write(7,*) 'Phase label in pix < 1--error at ',m
end if
if(pix(m).gt.nphase) then
write(7,*) 'Phase label in pix > nphase--error at ',m
end if
500 continue
return
end
c Subroutine to compute the total current in the x, y, and z directions
subroutine current(nx2,ny2,nz2,ns2,currx,curry,currz,u,gx,gy,gz)
complex u(ns2),gx(ns2),gy(ns2),gz(ns2),currx,curry,currz
complex cur1,cur2,cur3
c auxiliary variables involving the system size
nx=nx2-2
ny=ny2-2
nz=nz2-2
L22=nx2*ny2
c initialize the volume averaged currents
currx=cmplx(0.0,0.0)
curry=cmplx(0.0,0.0)
currz=cmplx(0.0,0.0)
c Only loop over real sites and bonds in order to get true total current
do 10 k=2,nz2-1
do 10 j=2,ny2-1
do 10 i=2,nx2-1
m=L22*(k-1)+nx2*(j-1)+i
c cur1, cur2, and cur3 are the currents in one pixel
cur1=0.5*( ( u(m-1)-u(m) )*gx(m-1)+( u(m)-u(m+1) )*gx(m) )
cur2=0.5*( ( u(m-nx2)-u(m) )*gy(m-nx2)+( u(m)-u(m+nx2) )*gy(m) )
cur3=0.5*( ( u(m-L22)-u(m) )*gz(m-L22)+( u(m)-u(m+L22) )*gz(m) )
c sum pixel currents into volume averages
currx=currx+cur1
curry=curry+cur2
currz=currz+cur3
10 continue
currx=currx/float(nx*ny*nz)
curry=curry/float(nx*ny*nz)
currz=currz/float(nx*ny*nz)
return
end