Create Straight Dislocations for ParaDiS Input
Jump to navigation
Jump to search
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