Create dislocations for ParaDiS input file
Sylvie Aubry (10/29/07)
implicit none integer nmax parameter (nmax = 20) integer ndis,n integer i,j,t,k,m,inc,m0,m1,method integer tmin(nmax),tmax(nmax),tstep(nmax) integer length,fp,maxstep integer ncount,ncount_int(nmax+1),ncount_tmp(nmax+1) integer seed(nmax) double precision bx(nmax),by(nmax),bz(nmax) double precision nx(nmax),ny(nmax),nz(nmax) double precision a(nmax),b(nmax),c(nmax) double precision x0(nmax),y0(nmax),z0(nmax) double precision x,y,z,tmp double precision rc,rTol double precision xmin,ymin,zmin,xmax,ymax,zmax double precision bnorm,nnorm double precision burgMag,L,Lmin,Lmax,minSeg,maxSeg character*50 name,name2,type
Read in some parameters
L : Length in micrometer
read(*,*) L
maxstep for the simulation
read(*,*) maxstep
L = (read-in L) micrometer = 10^-6 m.
burgMag = 2.725d-10 L = L*1d-6/burgMag Lmin = -L/2.d0 + 300.d0 Lmax = L/2.d0 write(*,*) write(*,*) 'L=',L write(*,*) 'Lmin=',Lmin write(*,*) 'Lmax=',Lmax write(*,*) 'Total number of steps=',maxstep
number of points on each dislocation line
inc = 50 write(*,5) 'Number of points on the curve=',inc write(*,*)
File name
read(*,*) name
ndis = how many dislocations
read(*,*) ndis if (ndis.le.0.and.ndis.gt.nmax) then write(*,*) 'Wrong number of dislocations.' stop else write(*,*) 'Reading ',ndis,' dislocations.' endif
Read Burgers vectors and glide planes
do i=1,ndis write(*,*) write(*,30) 'For dislocation ',i,'/',ndis
Burgers vector
read(*,*) bx(i),by(i),bz(i)
Glide planes
read(*,*) nx(i),ny(i),nz(i)
Normalize Burgers vector/Glide plane in units of burgMag = a sqrt(3)/2.
do j=1,3 bnorm = sqrt(bx(i)**2 + by(i)**2 + bz(i)**2) bx(i)= bx(i)/bnorm by(i)= by(i)/bnorm bz(i)= bz(i)/bnorm
nnorm = sqrt(nx(i)**2 + ny(i)**2 + nz(i)**2) nx(i)= nx(i)/nnorm ny(i)= ny(i)/nnorm nz(i)= nz(i)/nnorm enddo
Choose how to get dislocation line
read(*,*) method if (method.eq.0) then
Define dislocation line automatically according to type and seed
Type of the dislocation : screw, edge, mixed
read(*,*) type write(*,*) 'Type of dislocation: ',type(1:length(type)) read(*,*) seed(i) call get_abc(a(i),b(i),c(i),x0(i),y0(i),z0(i), & bx(i),by(i),bz(i),nx(i),ny(i),nz(i), & type,Lmin,Lmax,seed(i)) else read(*,*) a(i),b(i),c(i),x0(i),y0(i),z0(i) endif
write(*,20) 'bx=',bx(i),' by=',by(i),' bz=',bz(i) write(*,20) 'nx=',nx(i),' ny=',ny(i),' c=',nz(i)
write(*,20) 'a=',a(i),' b=',b(i),' c=',c(i) write(*,20) 'x0=',x0(i),' y0=',y0(i),' z0=',z0(i)
call get_bounds(xmin,xmax,Lmin,Lmax,a(i),x0(i)) call get_bounds(ymin,ymax,Lmin,Lmax,b(i),y0(i)) call get_bounds(zmin,zmax,Lmin,Lmax,c(i),z0(i))
tmin(i) = max(xmin,ymin,zmin) tmax(i) = min(xmax,ymax,zmax)
if (tmin(i).gt.tmax(i)) then tmp = tmin(i) tmin(i) = tmax(i) tmax(i) = tmp endif
inc in the final number of points on the dislocation line.
tstep(i) = (tmax(i)-tmin(i))/inc if (tstep(i).eq.0) then write(*,*) 'Found a number of point on the ', & 'curve equals to zero - ABORT' stop endif write(*,*) 'bounds=',tmin(i),tmax(i) enddo
Write outputs files
write .cn file
name2 = 'Outputs/'//name(1:length(name))//'_results' open(unit=20,file='../../Runs/'//name(1:length(name))//'.cn') write(20,400) '### ParaDis input file (create_dislocation.f) ###' write(20,400) '### S.A. (10/29/07)' write(20,*) write(20,400) '#Directory to write output files' write(20,*) 'dirname =',name2(1:length(name2)) write(20,*) write(20,400) '#Input files for PBC image stresses' write(20,400) 'Rijmfile = "Inputs/Rijm.cube.out"' write(20,400) 'RijmPBCfile = "Inputs/RijmPBC.cube.out"' write(20,*) write(20,400) '#Total simulation step' write(20,*) 'maxstep=',maxstep write(20,*) write(20,400) '#The total number of CPUs' write(20,400) 'numXdoms = 1' write(20,400) 'numYdoms = 1' write(20,400) 'numZdoms = 1' write(20,*) write(20,400) '#Cells for dislocation grouping (>=3)' write(20,400) 'numXcells = 3' write(20,400) 'numYcells = 3' write(20,400) 'numZcells = 3' write(20,*)
write(20,400) '#Fast multipole method specs.' write(20,*) 'fmEnabled = 0 #disable fast multipole' write(20,*) 'fmMPOrder = 2' write(20,*) 'fmTaylorOrder = 4' write(20,*) 'fmCorrectionTbl = "Inputs/fm-ctab.m2.t4.dat"' write(20,*)
write(20,400) 'timestepIntegrator = "backward-euler"' write(20,*) write(20,400) '#Simulation box size' write(20,*) 'xBoundMin =',Lmin write(20,*) 'xBoundMax =',Lmax write(20,*) 'yBoundMin =',Lmin write(20,*) 'yBoundMax =',Lmax write(20,*) 'zBoundMin =',Lmin write(20,*) 'zBoundMax =',Lmax write(20,*)
write(20,400) '#Boundary conditions' write(20,*) 'xBoundType = 0' write(20,*) 'yBoundType = 0' write(20,*) 'zBoundType = 0' write(20,*)
write(20,400) '#Mobility law function' write(20,*) 'mobilityLaw = "BCC_0"' write(20,*) 'MobScrew =', 10. write(20,*) 'MobEdge =', 10. write(20,*) 'MobClimb =', 1.e-4
max/min Seg = max/min length of dislocation segment. current dislocation segment length is inc
write(20,*) minSeg = 4*inc maxSeg = 2*minSeg write(20,400) '#Discretization' write(20,*) 'maxSeg =',maxSeg write(20,*) 'minSeg =',minSeg write(20,*)
write(20,400) '#Maximum nodal displacement at each time step' write(20,400) 'rmax = 100' write(20,*) write(20,400) '#Error tolerance in determining time step'
It should be such that rntol < 0.5 rc = 0.5 a.
rc = minSeg/sqrt(3.d0)*0.5 rTol = 0.45*0.5*rc write(20,*) 'rTol =',rTol write(20,*) write(20,400) '#Maximum time step' write(20,*) 'maxDT = 1e+07' write(20,*) write(20,*) '#Core cut-off radius' write(20,*) 'rc =',rc write(20,*)
write(20,400) '#Turn on elastic interaction' write(20,*) 'elasticinteraction = 1' write(20,*)
write(20,400) '#Applied stress in Pa (xx,yy,zz,yz,zx,xy)' write(20,*) 'appliedStress = [ 0 0 0 0 0 0 ]' write(20,*)
write(20,400) '#Loading' write(20,*) 'loadType = 1' write(20,*) 'eRate = 1' write(20,*) 'edotdir = [1 0 0]' write(20,*)
write(20,400) '#Save files' write(20,*) 'savecn = 1' write(20,*) 'savecnfreq = 100' write(20,*)
write(20,400) '#Povray output' write(20,*) 'povray = 1' write(20,*) 'povrayfreq = 100' write(20,*) 'povraydt = -1.000000e+00' write(20,*) 'povraytime = 0.000000e+00' write(20,*) 'povraycounter = 0' write(20,*) 'saveprop = 1' write(20,*) 'savepropfreq = 1' write(20,*) write(20,*) 'preserveOldTags = 1' close(20)
Write .data file
open(unit=20,file='../../Runs/'//name(1:length(name))//'.data') write(20,400) '### ParaDis nodal file (create_dislocation.f)' write(20,400) '### S.A. (10/29/07)###' write(20,400)
write(20,400) '#File version number' write(20,400) '2' write(20,400)
write(20,400) '#Number of data file segments' write(20,400) '1' write(20,400) '#Minimum coordinate values (x, y, z)' write(20,*) -L/2,-L/2,-L/2 write(20,400) '#Maximum coordinate values (x, y, z)' write(20,*) L/2,L/2,L/2 write(20,*)
write(20,400) '#Node count' ncount = 0 do i=1,ndis+1 ncount_int(i) = 0 enddo
do i=1,ndis do t=tmin(i),tmax(i),tstep(i) ncount = ncount + 1 ncount_int(i+1) = ncount_int(i+1) + 1 enddo enddo
write(20,*) ncount write(20,400) '#Domain geometry (x, y, z)' write(20,400) '1 1 1' write(20,*) write(20,400) '#Domain boundaries (x, y, z)' write(20,*) -L/2 write(20,*) -L/2 write(20,*) -L/2 write(20,*) L/2 write(20,*) L/2 write(20,*) L/2 write(20,*)
write(20,400) '#node_tag, x, y, z, num_arms, constraint' write(20,400) '#arm_tag, burgx, burgy, burgz, nx, ny, nz' write(20,400) '#length in unit of burgMag' m = 0 m0 = 0 m1 = 0 do i=1,ndis
first
m0 = m m1 = m1 + ncount_int(i) t=tmin(i) x = x0(i) + a(i)*t y = y0(i) + b(i)*t z = z0(i) + c(i)*t write(20,250) '0,',m,' ',x,y,z,2,0 write(20,300) ' 0,',m0+ncount_int(i+1)-1, & bx(i),by(i),bz(i) write(20,350) ' ',nx(i),ny(i),nz(i) write(20,300) ' 0,',m+1,-bx(i),-by(i),-bz(i) write(20,350) ' ',nx(i),ny(i),nz(i) m = m + 1
all
do t=tmin(i)+tstep(i),tmax(i)-tstep(i),tstep(i) x = x0(i) + a(i)*t y = y0(i) + b(i)*t z = z0(i) + c(i)*t write(20,250) '0,',m,' ',x,y,z,2,0 write(20,300) ' 0,',m-1,bx(i),by(i),bz(i) write(20,350) ' ',nx(i),ny(i),nz(i) write(20,300) ' 0,',m+1,-bx(i),-by(i),-bz(i) write(20,350) ' ',nx(i),ny(i),nz(i) m = m + 1 enddo
last node
t = tmax(i) x = x0(i) + a(i)*t y = y0(i) + b(i)*t z = z0(i) + c(i)*t write(20,250) '0,',m,' ',x,y,z,2,0 write(20,300) ' 0,',m-1,bx(i),by(i),bz(i) write(20,350) ' ',nx(i),ny(i),nz(i) write(20,300) ' 0,',m1,-bx(i),-by(i),-bz(i) write(20,350) ' ',nx(i),ny(i),nz(i) m = m + 1 enddo
close(20)
5 format(a,i3) 10 format(a,f9.2,a) 20 format(a,f9.2,a,f9.2,a,f9.2) 30 format(a,i2,a,i3) 100 format(3f15.3,i2,a) 200 format(2i5,6f10.3,a)
250 format(a,i4,a,3f12.2,2i2) 300 format(a,i4,3f12.2) 350 format(a,3f12.2) 400 format(a)
end
double precision function min(a,b,c) implicit none double precision a,b,c,min0 if (a.lt.b) then min0 = a else min0 = b endif if (min0.lt.c) then min = min0 else min = c endif return min end
double precision function max(a,b,c) implicit none double precision a,b,c,max0 if (a.gt.b) then max0 = a else max0 = b endif if (max0.gt.c) then max = max0 else max = c endif end
integer function length(str) character*(*) str n = len(str) 10 if (n.gt.0) then if (str(n:n).eq.' ') then n = n - 1 goto 10 endif endif length = n return end
subroutine random_gen(ranval,min,max,seed) integer seed real*4 ranval_tmp double precision min,max,ranval ranval_tmp = rand(seed) ranval = ranval_tmp*(max-min) + min end
subroutine get_abc(a,b,c,x0,y0,z0,bx,by,bz,nx,ny,nz, & type,Lmin,Lmax,seed) implicit none integer length integer seed,seedX0,seedY0,seedZ0,seedalpha double precision a,b,c double precision x0,y0,z0,bx,by,bz,nx,ny,nz,Lmin,Lmax double precision alpha character*20 type
seedX0 = seed + 19 seedY0 = seed + 395 seedZ0 = seed + 34 seedalpha = seed
call random_gen(x0,Lmin/5d0,Lmax/5d0,seedX0) call random_gen(y0,Lmin/5d0,Lmax/5d0,seedY0) call random_gen(z0,Lmin/5d0,Lmax/5d0,seedZ0) call random_gen(alpha,-100.d0,100.d0,seedalpha)
screw : dislocation line is parallel to the Burgers vector.
if (type(1:length(type)).eq.'screw') then a = alpha * bx b = alpha * by c = alpha * bz endif
edge : dislocation line is normal to the Burgers vector and the glide plane
if (type(1:length(type)).eq.'edge') then a = alpha * (by*nz - bz*ny) b = -alpha * (bx*nz - bz*nx) c = alpha * (bx*ny - by*nx) endif
mixed : edge and screw component
if (type(1:length(type)).eq.'mixed') then a = alpha * bx + alpha * (by*nz - bz*ny) b = alpha * by - alpha * (bx*nz - bz*nx) c = alpha * bz + alpha * (bx*ny - by*nx) endif
end
subroutine get_bounds(xmin,xmax,Lmin,Lmax,a,x) implicit none integer i double precision xmin,xmax,Lmin,Lmax,a,x if (a.ne.0) then if (a.gt.0) then xmin = (Lmin-x)/a xmax = (Lmax-x)/a else xmin = (Lmax-x)/a xmax = (Lmin-x)/a endif else xmin = -10d10 xmax = 10d10 endif end
subroutine get_rand(seed,yy) implicit none integer i, number, old, seed, x, y double precision yy number = 1000 old = seed do i = 1, number x = MOD((57*old+1), 256) y = MOD((57*x+1), 256) old=y enddo yy = dabs(y*1.d0/255.d0) end