DDLab Manual 01 - Revision history

Caiwei at 18:01, 21 December 2007

2007-12-21T18:01:27Z

← Older revision		Revision as of 18:01, 21 December 2007
Line 22:		Line 22:
	* Time controls : totalstpes, dt0		* Time controls : totalstpes, dt0
	* Display controls : plotfreq, plim, viewangle, printfreq, printnode		* Display controls : plotfreq, plim, viewangle, printfreq, printnode
	* Materials constant : MU, NU, Ec, a		* Materials constant : MU (shear modulus), NU (Poisson's ratio), Ec (core energy), a (core radius)
	* Miscellaneous : appliedstress, maxconnections		* Miscellaneous : appliedstress, maxconnections

Caiwei: mex compiler

2007-12-19T19:09:50Z

mex compiler

New page

== Introduction ==
The purpose of this document is to give a brief description of [http://micro.stanford.edu/~caiwei/Forum/2005-12-05-DDLab/ DDLab]. If advanced users are interested in the algorithms of
the code, they may simply look at the source code. We have tried to write and organize
the subroutines in a compact yet readable manner.

<font color="red">Warning</font>: The main executable file <tt>dd3d.m</tt> requires the compilation of

mex -O SegSegForces.c

You need to have a C compiler installed on your computer to make it work. This is ok when you use <tt>gcc-3.3</tt>. You may have difficulty compiling with higher version of gcc.

== Input parameters of DDLab ==
There are various input parameters in DDLab. In this section, we don't give a detailed expla-
nation of each parameter but we explain some. The relation between parameters will
be explained in more detail in the following sections.

We can divide all input parameters into the following categories.

* Dislocation structure : rn, links* Mobility : mobility
* Integrator : integrator
* Topological changes : lmax, lmin, areamin, areamax, rmax, rann, rntol, doremesh, docollision, desepraration
* Time controls : totalstpes, dt0
* Display controls : plotfreq, plim, viewangle, printfreq, printnode
* Materials constant : MU, NU, Ec, a
* Miscellaneous : appliedstress, maxconnections

== Nodal representation of dislocation structure==
Figure 1: Nodal representation of dislocation structure
Fig. 1 shows a simple approach that can represent an arbitrary dislocation network. The
dislocations are specified by a set of nodes that are connected with each other by straight
segments. Each segment has a nonzero Burgers vector. Because the Burgers vector is
defined only after a sense of direction is chosen for the dislocation line, we can define bij

== Mobility ==
In DDLab, the motion of dislocation is expressed by the motion of nodes. Thus,
mobility gives the velocity. DDLab contains several mobility
laws; mobfcc1.m, mobfcc0.m, and mobbcc0.m.

==Integrator==
An integrator scheme is used to calculate the position of each nodes for a given time step. The nodal
forces and mobility model comprise the equation motion for the nodes
<math> v_i = \frac{dr_i}{dt} = g_i(r_j)g</math>
where g_i sums up both the nodal force and the mobility model. This is a
first ordinary differential equation and we solve it by integration. DDLab contains two numerical integrator:
<tt>eulerforward.m</tt>, <tt>newtonkrylov.m</tt>. <tt>ode15s.m</tt>,
<tt>trapezoid.m</tt>. In this work, we use <tt>trapezoid.m</tt> is used.

==Topological Changes==
In order to accurately represent dislocation behavir and also for numerical efficiency,
we need a way to handle topological changes, i.e.
changes in the connectivity of nodes. For example, should a line change its length and/or
curvature during the simulation, it may become necessary to adjust the number of nodes
used to represent this dislocation line. The following parameters are used for topological
changes in DDLAB.
* lmax: Maximum length of a dislocation segment
* lmin: Minimum length of a dislocation segment
* areamin: Maximum area criterion
* areamax: Minimum area criterion
* rmax: maximum distance a node may travel in one cycle
* rann: annihilation distance used to calculate the collision of dislocation lines
* rntol: solution tolerance used to control the automatic timestepping
* doremesh: Toggle remesh.m for remeshing nodes. doremesh is set as zero to disable
the use of remesh rule. When it is set as zero, it is enabled.

=Time controls=
The following parameters are for time controls.

* totalsteps : number of cycles that are run for completion of dd3d command
* dt0 : largest timestep that can be taken during a cycle

=Display controls=
The following parameters are for display controls

* plotfreq : number of cycles between monitored node write statements
* plim : limits of plotting space
* viewangle : angle of viewpoint for the 3D plot of the geometry
* printfreq : number of cycles between monitored node write statements
* printnode : nodeid of the monitored node

=Materials constants=
The following parameters are Materials constants used in DDLAB.
* Mu: Shear modulus
* NU: Poisson's ratio
* Ec: Dislocation core energy
Dislocation core radius used for non-singular force calculation

=Miscellaneous=
* maxconnections : maximum number of segments a node may have

==Flow of the code in DDLab==
After setting up the initial configuration and other parameters, both DDLab and ParaDiS will follow the same
algorithm to simulate dislocation dynamics, as is listed below:

* calculate the force of each node
* calculate the velocity of each node
* calculate the new position of each node
* make necessary topological changes
* repeat (1) to (4) until the maximum step is reached

DDLab just track all nodal positions and plot them. However, to
better understand the code structure, two
flow charts are presented in this primer:
subroutine map (Fig. 2) and function map (Fig. 3). The subroutine map consists of the
MATLAB subroutine file name, but the function map consists of the function of each
subroutine file. Thus, if you want to know the code in detail, just compare the subroutine
map with the function map, then find the file you want and take a look at. Thus, these two
maps are very useful to study the code structure of DDLab. All functions of subroutines are
explained in the following section which includes some physics for the detailed explanations.

==Subroutines details==
* <tt>SegSegforce.c</tt>

This subroutine calculates the force on the two end nodes between two dislocation seg-
ments (so, four nodes) in Fig. 4. The sum of these four force vectors should be zero by
action-reaction principle.
<math>f1 + f2 + f3 + f4 = 0</math>

Figure 4: The interaction between the two dislocation segments

* <tt>cleanupnodes.m</tt>
This subroutine start by removing all of the bad nodes in the system and cleans up the
data structures to include good nodes only. A bad node is a point without connection
to any dislocation segment. The point is usually generated by the shrinkage of a dislocation
loop as shown in Fig. 5. At the final stage, the loop consists of two nodes, so the line
direction of two segments is opposite, but their Burgers vector is the same. Thus, two segments
should annihilate. Then, only two nodes remains. These two nodes are
only a mathematical points in elastic medium. In real material, the loop disappears after
shrinkage, so the two nodes in DDLab has to be removed. Thus, <tt>cleanupnodes.m</tt> performs
the removal of these kinds of isolated nodes.

Note:

This is very similar with to another case (two nodes and two segments) of <tt>collision.m</tt>.
Sometimes segments are not permitted to annihilate when the sum of Burgers vectors are not zero.
However, in the case of loop shrinkage, two segments always have the opposite direction of
line sense with the same Burgers vector. Thus, the annihilation occurs. See <tt>collision.m</tt>
section.

* <tt>genconnectivity.m</tt>
This subroutine generates connectivity and <tt>linksinconnect</tt>. In fact, rn and links
matrices are enough to describe the dislocation structure. These two matrices give <tt>linkid</tt>
when we know <tt>nodeid</tt>. Conversely, <tt>connectivity</tt> and <tt>linksinconnect</tt> give <tt>nodeid</tt> when

Figure 5 The shrinkage of a dislocation loop.
we know linkid Here, nodeid is the node number which is the row number of rn matrix.
linkid is the link number which is the row number of links matrix.

==Integrator==
Integrator calls <tt>drndt.m</tt>, <tt>drndt.m</tt> computes velocity and force on each node, and
force on each segment. Then, by using these values, Integrator integrates the velocity
with respect to time step. Then, we can get the nodal positions at the current time step.
The nodal forces and the mobility model comprise the equations of motion of nodes:

<math>
v_i = \frac{dr_i}{dt} = g_i(r_j)g
</math>

where function <math>g_i</math> sums up both the nodal force and the mobility model. This is first
ordinary differential equation (ODE), for which the simplest numerical integrator is the so-
called Euler forward method:

<math>
ri(t + �t) = ri + gi(rj(t))�t:
</math>

A numerical integrator with a higher order of accuracy is the trapezoid method,
<math> r_i(t + �t) =
g_i(r_j(t)) + g_i(r_j(t + �t))
</math>

Another important feature of integrator is the time step adjustment. When the time
step is large, we can perform simulation which require long physical time. However, the
large time step always bring large error. Conversely, when the time step is too short, we can
perform the more precise simulation, but it take too long CPU time to get the result. Thus,
it is better to adjust the time step in every cycle. The present version of DDLab has the
time step adjustment as below. This is the part of the integrator code.

...
err=rnvec1-rnvec0-(vnvec+vnvec0)/2*dt;
...
errmag=max(abs(err));
...
if(errmag<rntol)
break;
else
dt=dt/2;
...
maxchange=1.2;
...

If the defined err is larger than rntol, the time step is adjusted as the half of current
time step. If the defined err is smaller than rntol, the time step is adjusted as 1.2 times of
current time step. Thus, DDLab changes its time step dynamically, enabling a much more efficient
simulation.

A series of time integration schemes are also included with the distribution named
<tt>int eulerforward</tt>, <tt>int ulerbackward</tt>, <tt>int newtonkrylov</tt>.
<tt>int eulerforward is an explicit forward integration scheme. <tt>int eulerbackward</tt> is a simple implicit time
integration scheme, and <tt>int newtonkylov</tt> is a more sophisticated inexact implicit time integration scheme based
on the Newton-Krylov matrix-free solution method. The core of the
<tt>int newtonkrylov</tt> employs a method written by C. T. Kelley. The time integration scheme
is specified in the input deck with the integrator input parameter.

4.0.14 drndt.m
drndt.m calls segforcevec.m. segforcevec.m calculates the force on each segment,
then return it to drndt.m. Then, drndt.m calls mobility. mobility calculates the force
and velocity on each node and return them to drndt.m Finally, drndt.m return the velocity
on node to integrator, which integrate the velocity with respect to time step. Then, we
can get the nodal positions at the current time step.
4.0.15 segforcevec.m
segforcevec calculates the force on each segment. Note that this is the force not on
the node but on the segment. It is because DDLAB use Peach-Koehler formula
fPK(x) = (�(x) � b) � �(x)
.
When we compute PK force, we need to know the line sense vector (�) for segment.
Thus, �rstly DDLAB calculates the force on each segment which has two end nodes. Then,
mobility function gives the half of force to each node. Finally, we can get the nodal
force. segforcevec.m includes the four functions; constructsegmentlist, pkforcevec,
selfforcevec, remoteforcevec. constructsegmentlist give us the list of all information
about segments as the format
n0, n1, bx, by, bz, x0, y0, z0, x1, y1, z1, nx, ny, nz,
where n0 and n1 are nodeids of two end nodes. bx, by, bz is the Burgers vector of one
segment which consists of two end node, n0 and n1. x0, y0, z0 is the position of node n0
and x1, y1, z1 is the position of node n1. nx, ny, nz is the normal vector of the glide
plane on which the dislocation segment lies.
Based on this list, segforcevec calculate the forces on each segment. selfforcevec
calculates the force on each segment by itself. remoteforcevec calculates the force on the
segment on each segment by the other segments.
Finally, the sum of these forces is the force on segment as below.
fseg=[fpk, fpk]*0.5+[fs0, fs1]+[fr0, fr1];
The �rst component of fseg (fpk, fs0, fr0) are the force vectors on node n0 and the
second ones are the force vectors on the node n1. The node n0 and n1 are the end nodes of
each segment as shown in Fig. 5.

Figure 7: (a) The example dislocation structure. (b) DDLAB calcuate the force on
the segment by segforcevec.m (c) Then, the half of the force is assigned to two
end nodes of each segment (d) The force on the node has several arms is just the
summation of forces on corresponding arms.
4.0.16 mobfcc1.m, mobfcc0.m, mobfcc1.m
A dislocation's response to forces may also depends on the line orientation and varies
signi�cantly from one material to another. Similar to the core energy, how dislocation
respond to forces is a function of the non-linear atomistic interactions in the dislocation
core and, as such, is beyond the scope of linear elasticity theory. The response functions
have to be imported into the line DD model in the form of external materials input, either
from atomistic simulations or from experiments
A series of mobility functions are included in the distribution named mobbcc0 mobbcc0b
and mobfcc1. mobbcc0 and mobbcc0b are intended to simulate bcc behavior in which
the the screw dislocations are able to glide in any plane normal to their Burgers vector.
The di�erence between them is slight and appears in the mobility of dislocations of mixed
character. mobfcc1 is intended to simulate fcc behavior in which the screw dislocation are
not able to cross-slip. The mobility function is chosen in the input deck with the mobility
input parameter.

4.0.17 separation.m
For numerical and physical reasons, line-DD simulations need to handle topological
changes, i.e. changes on the connectivity between nodes since we may want to adjust the
number of nodes that represent a dislocation line if the line gets longer or shorter during the
simulation, or when two dislocation lines meet in space, they may either annihilate or zip
together to form a junction, which also results in a change of nodal topology. Thus many
types of topological changes can be encountered in a line-DD simulation. Fortunately, since
we use a nodal representation here, all topological changes can be implemented through two
basic operators: split (one node split into two ) and merge (two nodes merge into one) (See
next section) and . The implementation of these two operators is straightforward - all one
needs to do is to make sure that at the end of the operation the Burgers vector sum rule
at every node and segment is still satis�ed, moreover, two nodes are either disconnected or
connected only once, and each node is connected with at least two other nodes, if a node
has no segment, it will be deleted.
For example, for a 4-arm node such as P in there are 3 di�erent ways to partition its
arms: (12)(34), (13)(24) and (14)(23). It is reasonable to expect that the way nature would
choose should be the one that gives rise to the maximum energy dissipation rate, which is
de�ned below.
Suppose an n-arm mode i stays intact (not splitting) and it feels a force fi and will move
at velocity vi . Then the local energy dissipation rate is,
_W
i = fi � vi
. Now suppose that node i splits into two nodes P and Q, such that node P retains 1, ...
,s of the original neighbors, and node Q retains the remaining neighbors. Let fP and fQ
be the forces on the two nodes and vP and vQ be their velocities given by the mobility
function. Then the local energy dissipation rate is,
WP_ Q = fP � vP + fQ � vQ
. If _ WPQ > _Wi, then node i prefers to split into two nodes P and Q instead of moving
as a single node. The energy dissipation rate can be computed for all possible (topological
distince) modes to split i. The mode with the highest energy dissipation rate is preferred.
If a node will split in next step, the two new nodes actually stay at the same location
as the parent node at the current step. Because the velocities of the two nodes is di�erent
(otherwise the node should not split), the two nodes will be move away from each other in
the next time step.

* collision.m

* remesh.m

==Self consistency in DDLab units==

The unit system of DDLab input file should be self-consistent. We recommend to input all variables in S.I. units.

For instance, stresses in Pascal, distances in meter, times in seconds and velocities in m/s.

The mobility law is such that
<math>
v = M (\sigma \cdot b)
</math>

with b in meter, <math>\sigma</math> in Pascal, and v in m/s. That makes the mobility parameter M in <math> Pa^{-1} s^{-1}</math>.

We often use the drag coefficient <math>B = M^{-1}</math> in the mobility law. The unit of <math>B</math> is <math> Pa\, s</math>

As another option, we can change the unit for distance (including both nodal positions and Burgers vectors) to nm, whereas the unit for stress is still Pa, the unit of time is still second, and the unit for mobility is still <math> Pa^{-1} s^{-1}</math>.

Caiwei: /* Self consistency in DDLab units */

2007-12-18T22:18:33Z

Self consistency in DDLab units

← Older revision		Revision as of 22:18, 18 December 2007
Line 68:		Line 68:
	with b in meter, <math>\sigma</math> in Pascal, and v in m/s. That makes the mobility parameter M in <math> Pa^{-1} s^{-1}</math>.		with b in meter, <math>\sigma</math> in Pascal, and v in m/s. That makes the mobility parameter M in <math> Pa^{-1} s^{-1}</math>.

	We often use the drag coefficient <math>B = M^{-1}</math> in the mobility law. The unit of <math>B</math> is <math> Pa s</math>		We often use the drag coefficient <math>B = M^{-1}</math> in the mobility law. The unit of <math>B</math> is <math> Pa\, s</math>

	As another option, we can change the unit for distance (including both nodal positions and Burgers vectors) to nm, whereas the unit for stress is still Pa, the unit of time is still second, and the unit for mobility is still <math> Pa^{-1} s^{-1}</math>.		As another option, we can change the unit for distance (including both nodal positions and Burgers vectors) to nm, whereas the unit for stress is still Pa, the unit of time is still second, and the unit for mobility is still <math> Pa^{-1} s^{-1}</math>.

Caiwei: /* Self consistency in DDLab units */

2007-12-18T22:17:55Z

Self consistency in DDLab units

← Older revision		Revision as of 22:17, 18 December 2007
Line 57:		Line 57:
	==Self consistency in DDLab units==		==Self consistency in DDLab units==

	~~Any~~ DDLab input file is self-consistent. It ~~is recommended~~ to input all variables in ~~the~~ S.I. units.		The unit system of DDLab input file should be self-consistent. We recommend to input all variables in S.I. units.

	For instance, stresses in Pascal, distances in meter, times in seconds and velocities in m/s.		For instance, stresses in Pascal, distances in meter, times in seconds and velocities in m/s.
Line 66:		Line 66:
	</math>		</math>

	with b in meter, <math>\sigma</math> in Pascal, and v in m/s. That makes M in <math> Pa^{-1} s^{-1}</math>.		with b in meter, <math>\sigma</math> in Pascal, and v in m/s. That makes the mobility parameter M in <math> Pa^{-1} s^{-1}</math>.

			We often use the drag coefficient <math>B = M^{-1}</math> in the mobility law. The unit of <math>B</math> is <math> Pa s</math>
	We can see from this little exercise that if instead of using meters as the unit for distances, we use nanometers instead. The mobility law does not change units.

			As another option, we can change the unit for distance (including both nodal positions and Burgers vectors) to nm, whereas the unit for stress is still Pa, the unit of time is still second, and the unit for mobility is still <math> Pa^{-1} s^{-1}</math>.

Caiwei: /* Introduction */

2007-12-18T22:09:45Z

Introduction

New page

== Introduction ==
Professor Wei Cai in the Department of Mechanical Engineering at Stanford
University with the team at the Lawrence Livermore National Laboratory has developed
the new dislocation simulation code, ParaDiS (Parallel Dislocation Simulator) and DDLAB
(Dislocation Dynamics Lab.). Dislocation simulation enables us to control of the simulation
conditions and to access the positions of all dislocation lines at each given time. Since
the time and length scales in dislocation simulation codes are on the order of seconds
and microns, respectively, which are similar with experiments, they
can show realistic insights that help answer fundamental questions in small-scale plasticity, such as motion and interaction of dislocations in small volumes. With the use of
massively parallel computers, ParaDiS can deal with the large number of dislocations in an
effcient way. However, the structure of cods of ParaDiS written by C is quite complicated,
code modifications are not easy so DDLAB, a simpler matlab version of ParaDiS was designed. Since
MATLAB has built-in commands for computation, DDLab has a much simpler
structure than that of ParaDiS.

Although DDLAB does not support the parallel computation, it is enough to deal with simple dislocation problems, for instance, consisting of
two or three dislocations. It was created to test new codes before implementation into ParaDiS.

The purpose of this document is to give a brief description of the DDLAB software
for end users who are interested in writing input files and running simulations using the
software. This document does not provide an exhaustive description of the algorithms
in the code, but rather gives a description of necessary elements of the input deck for the
code and the execution procedures. If advanced users are interested in the algorithms of
the code, they may simply look at the source code. We have tried to write and organize
the subroutines in a compact yet readable manner.

== Input parameters of DDLab ==
There are various input parameters in DDLab. In this section, we don't give a detailed expla-
nation of each parameter but we explain some. The relation between parameters will
be explained in more detail in the following sections.

We can divide all input parameters into the following categories.
1. Dislocation structure : rn, links
2. Mobility : mobility
3. Integrator : integrator
4. Topological changes : lmax, lmin, areamin, areamax, rmax, rann, rntol,
doremesh, docollision, desepraration
5. Time controls : totalstpes, dt0
6. Display controls : plotfreq, plim, viewangle, printfreq, printnode
7. Materials constant : MU, NU, Ec, a
8. Miscellaneous : appliedstress, maxconnections

== Nodal representation of dislocation structure==
Before seeing rn and links, we need to know how to represent dislocation structure in
DDLAB.
Figure 1: Nodal representation of dislocation structure
Fig 1 shows a simple approach that can represent an arbitrary dislocation network. The
dislocations are specified by a set of nodes that are connected with each other by straight
segments. Each segment has a nonzero Burgers vector. Because the Burgers vector is
defined only after a sense of direction is chosen for the dislocation line, we can define bij

==Self consistency in DDLab units==

Any DDLab input file is self-consistent. It is recommended to input all variables in the S.I. units.

For instance, stresses in Pascal, distances in meter, times in seconds and velocities in m/s.

The mobility law is such that
<math>
v = M (\sigma \cdot b)
</math>

with b in meter, <math>\sigma</math> in Pascal, and v in m/s. That makes M in <math> Pa^{-1} s^{-1}</math>.

We can see from this little exercise that if instead of using meters as the unit for distances, we use nanometers instead. The mobility law does not change units.