In [ ]:
# import anything we need in this notebook
using MsUtils, Compose, Colors, PGFPlots
import LujiaLt; FEM = LujiaLt.FEM
In this notebook we collect the computational examples associated with Chapters 3 and 4 of Multiscale Simulation.
Remarks:
The codes are reasonably well structured, but not as far as a long-lived code would need to be. Rather, they are intended to quickly get some tests to run. See the main LujiaLt.jl library for a more sophisticated set of codes but still focusing in mathematical toy-models and Atoms.jl for the real thing.
These codes use some elementary functional style constructs.
These may initially require some getting used to, but in the long run they
make for code that is brief, readable and less likely to contain bugs.
In particular, Julia encourages to avoid if statements and instead use
its dispatch mechanism. This is a paradigm that makes for simple and
elegant code.
Compose.jl; to produce PDFs the latest master needs to be checked outColors.jlPGFPlots.jl: for beautiful Latex look-and-feel graphics.The homogeneous lattice model is $$ E(u) = \sum_{\ell \in \Lambda^{\rm hom}} \sum_{i = 1}^3 \phi\big( D_{a_i} u(\ell) \big) = \sum_{b \in B^{\rm hom}} \phi( Du_b ). $$
We can introduce a point defect by introducing additional lattice sites and associated bonds. A screw dislocation is introduced by adding a pre-strain to the bonds. In order to implement only one model, we will describe both of these within the same framework.
Let $\Lambda$ be a point defect reference configuration and $B$ an associated set of nearest-neighbour bonds. Specifically we consider two cases:
The pre-strain is a map ${\sf e} : B \to \mathbb{R}$. The energy is now given by $$ E(u) = \sum_{b \in B} \phi( {\sf e}_b + Du_b ) - \phi( {\sf e}_b ). $$
We define the two benchmark problems through specifying $\Lambda, B, {\sf e}$:
More details about these choices can be found in Chapters 3 and 4 of the book.
We first consider the purely atomistic approximation where only the admissible displacements are restricted. For $R \in \mathbb{N}$ and a defect core position $\xi$ the computational domain is given by $$ \Omega_R = \big( \Lambda \cap B_R(\xi) \big) - \xi $$ with Dirichlet conditions (clamped) on the boundary sites.
In [ ]:
# [1] Preliminaries
# a nice Julia hack that allows us to use dispatch effectively
NoDefectT = Val{:none}; const NoDefect = NoDefectT()
ScrewT = Val{:screw}; const Screw = ScrewT()
InterstitialT = Val{:int}; const Interstitial = InterstitialT()
NoneOrInt = Union{NoDefectT, InterstitialT}
NoneOrScrew = Union{NoDefectT, ScrewT}
"distance of several points from a centre"
dist(X::Matrix, x0) = dist(X .- x0)
dist(X::Matrix) = sqrt(sumabs2(X,1))
# default interstitial position
const ξi = [0.5; 0.0]
const ξs = [0.5; sqrt(3)/6]
# off-centre defect positions, to break symmetry
const ξi_off = [2/5; 1/10]
const ξs_off = [1/4; 1/8]
"default defect cores"
defcore(::InterstitialT) = ξi
defcore(::ScrewT) = ξs
defcore(::NoDefectT) = [0.0;0.0]
# [1] computional domain ------------------------------------------
"""computational domain; returns (x, y) coordinate
vectors of lattice points"""
domain(R, defect::NoneOrScrew; ξ=defcore(defect)) =
LujiaLt.lattice_ball(R=R, A=LujiaLt.Atri, x0=ξ)[1] .- ξ
domain(R, ::InterstitialT; ξ=ξi) = [domain(R, NoDefect, ξ=ξ) [0.0;0.0]]
# [3] bonds and operations on bonds -------------------------------
"If v : X → R, then fdiff(v) : Bonds → R storing Dv_b"
fdiff(v::Vector, B::Tuple) = v[B[2]] - v[B[1]]
fdiff(v::Matrix, B::Tuple) = v[:, B[2]] - v[:, B[1]]
fdiff(v, B::Matrix) = fdiff(v, mat2tup(B))
"""
bonds within a computational domain: find all the Delaunay
edges with length <= 1, or connected to the defect site.
"""
function bonds(X)
E = FEM.edges(FEM.Triangulation(X))
I = find( (sumabs2(fdiff(X, E), 1) .<= 1.0+1e-10)[:] .+
(min(dist(X[:, E[1,:][:]]), dist(X[:, E[2,:][:]]))[:] .<= 1e-10) )
return mat2tup(E[:, I])
end
"bond midpoints"
midpoints(X, B) = 0.5*(X[:, B[1]] + X[:, B[2]])
"distance of bond mid-points from x0"
bonddist(X::Matrix, B::Tuple, ξ) = dist( midpoints(X, B), ξ)
bonddist(X::Matrix, B::Tuple) = bonddist(X, B, [0.0; 0.0])
# [4] prestrain is a function of lattice positions ----------------------------
# for the interstitial case, instead of e_b = 1 - len(b)
# we take a signed version of this to ensure that the
# choice of direction of bond does not influence the result.
"compute prestrains"
prestrain(X, B, ::NoneOrInt) =
((1.0 - sqrt( sumabs2(fdiff(X, B), 1) )) .* sum(fdiff(X, B),1))[:]
predictor(X, ::ScrewT) = arg(X[1,:][:], X[2,:][:]) / (2*π)
prestrain(X, B, ::ScrewT) = grad2strain(fdiff(predictor(X, Screw), B))
arg(x, y) = angle(x + im * y)
grad2strain(d) = mod(d + 0.5, 1) - 0.5
# [5] last aspect of geometry: free indices ----------------------------------
Ifree(X, R) = find(dist(X) .< R)
Iclamp(X, R) = find(dist(X) .>= R)
;
In [ ]:
# Interstitial
X = domain(4.0, Interstitial; ξ=ξi)
B = bonds(X)
ee = prestrain(X, B, Interstitial)
MsUtils.plot_int(X, B, ee, plotwidth=10cm);
In [ ]:
E = FEM.edges(FEM.Triangulation(X))
X[:, E[1,:][:]]
In [ ]:
# Screw: the red dot denotes the core position
X = domain(4.0, Screw, ξ=ξs)
B = bonds(X)
ee = prestrain(X, B, Screw)
MsUtils.plot_screw(X, B, ee, plotwidth=10cm);
In [ ]:
# [6] Model type : ------------------------------------------------
# to collect all the information
# the construction is produced below
type AtmModel{T}
R # domain parameter
X # lattice coordinates
B # bonds
e # pre-strains
Ifree # free dofs
defect::T # type of defect
end
function AtmModel(; R=5, defect = NoDefect, ξ=defcore(defect))
X = domain(R, defect, ξ=ξ)
B = bonds(X)
return AtmModel(R, X, B, prestrain(X, B, defect), Ifree(X, R-1.1), defect)
end
"number of lattice sites in the computational domain"
nsites(m) = size(m.X, 2)
# [7] Energy ------------------------------------------------------
# * we could use ForwardDiff.jl to implement the gradient and hessian
# automatically, but it is very slow it is not difficult to do it by hand.
# * there are several choices for the pair potential, we keep one and leave
# the others commented out
##### OPTION 1: basic periodic pair potential : problem is it has ϕ'''(0) = 0
# "NN pair potential (we use the same potential for all experiments)"
# ϕ(r) = sin(r*π).^2
# "first derivative of ϕ"
# Dϕ(r) = π * sin(2*r*π)
# "second derivative of ϕ"
# D²ϕ(r) = 2*π^2 * cos(2*r*π)
##### OPTION 2: a quartic polynomial; this is the most robust choice
const CC = (0.1, 1.0, 0.333, 5.678)
"NN pair potential (we use the same potential for all experiments)"
ϕ(r) = CC[1]*r + 0.5*CC[2]*r.^2 + (CC[3]/6)*r.^3 + (CC[4]/24)*r.^4
"first derivative of ϕ"
Dϕ(r) = CC[1] + CC[2]*r + (CC[3]/2)*r.^2 + (CC[4]/6)*r.^3
"second derivative of ϕ"
D²ϕ(r) = CC[2] + CC[3]*r + (CC[4]/2) * r.^2
# ##### OPTION 3: modified periodic pair potential with ϕ'''(0) ≠ 0
# ##### with this one, the solver doesn't converge
# const C3 = 0.1
# const C1 = 1.0
# "NN pair potential (we use the same potential for all experiments)"
# ϕ(r) = C1*sin(r*π).^2 + C3 * sin(2*r*π).^3
# "first derivative of ϕ"
# Dϕ(r) = π * C1*sin(2*r*π) + 3*C3*π*sin(2*π*r).*sin(4*π*r)
# "second derivative of ϕ"
# D²ϕ(r) = π^2 * ( C1*2.0 * cos(2*r*π) - C3*3*sin(2*π*r) + C3*9*sin(6*π*r) )
"potential energy (difference), as function of displacement"
energy(m::AtmModel, U) = sum( ϕ(m.e + fdiff(U[:], m.B)) - ϕ(m.e) )
"gradient of potential energy"
gradient(m::AtmModel, U) = _grad_( m.B, Dϕ(m.e + fdiff(U[:], m.B)), length(U) )
_grad_(B, dϕ, N) = binsum([dϕ; -dϕ], [B[2]; B[1]], N)
"hessian as function of displacement"
hessian(m::AtmModel, U) = _hess_( m.B, D²ϕ(m.e + fdiff(U[:], m.B)), length(U) )
_hess_(B, hϕ, N) = sparse( [B[1]; B[2]; B[1]; B[2]],
[B[1]; B[2]; B[2]; B[1]],
[ hϕ; hϕ; -hϕ; -hϕ], N, N )
# [8] Basic Newton solver --------------------------------------------------
# this is normally not robust enough, but sufficient for our purposes.
function solve(m; x = zeros(nsites(m)), show = false)
If = m.Ifree
show && @printf("------Newton Iteration------\n")
nit = 0
for nit = 0:15
∇E = gradient(m, x)
show && @printf("%d : %4.2e \n", nit, norm(∇E[If], Inf))
norm(∇E[If], Inf) < 1e-8 && break
x[If] = x[m.Ifree] - hessian(m, x)[If, If] \ ∇E[If]
end
nit == 12 && warn("the Newton iteration did not converge")
show && @printf("----------------------------\n")
return x
end
;
We perform some basic tests to make sure the implementation is correct. The basic philosophy of these tests is that, if the energy, gradient and hessian are all consistent then the implementation is likely correct. Hense we perform a finite difference test. We will test all partial derivatives for decreasing steps $h$.
Note that $f'(x) = \frac{f(x+h) - f(x)}{h} + O(h)$, but when we take floating-point errors into account, then we obtain $$ f'(x) = \frac{f(x+h) - f(x)}{h} + O\left( h + \frac{\epsilon}{h} \right), $$ where $\epsilon$ denotes floating point accuracy. Hence, as $h$ decreases we first see an improvement but eventually a deterioration of the test agreement.
These tests use fdtest, which is implemented in MsUtils.jl.
In [ ]:
# run the FD test : change `defect` to :none, :int, :screw
# to test the three implementations
m = AtmModel(R = 10, defect = Interstitial)
println("Test Energy as function of displacement")
fdtest(0.1 * rand(nsites(m)), x->energy(m,x), x->gradient(m,x), x->hessian(m,x))
Our second test concerns the Newton scheme. If it is correctly implemented and if gradient and hessian are consistent (we just tested this), then it should converge in 4-5 iterations at most.
In [ ]:
# change `defect` to {NoDefect, Interstitial, Screw} for the three model problems
m = AtmModel(R = 10, defect=Interstitial)
solve( m; show=true, x = 0.0 * rand(nsites(m)) );
In [ ]:
# interstitial
m = AtmModel(R = 20, defect=Interstitial)
u = solve(m)
Du = fdiff(u, m.B)
axis = [ξi[1]-7; ξi[1]+7; ξi[2]-5; ξi[2]+5]
MsUtils.plot_strain(m.X, m.B, log(abs(m.e+Du)+2e-4);
axis=axis, cmap = colormap("blues"));
In [ ]:
# screw dislocation
m = AtmModel(R = 20, defect=Screw)
u = solve(m)
Du = fdiff(u, m.B)
axis = [ξs[1]-7; ξs[1]+7; ξs[2]-5.7; ξs[2]+4.3]
MsUtils.plot_strain(m.X, m.B, log(abs(m.e+Du)+3e-2);
axis=axis, cmap = colormap("blues") )
# filename = "apx-strain-screw.svg", printwidth=12cm);
In [ ]:
# ENVELOPES of strain fields
R = 50 # computational domain radius
Rp = 20 # radius of domain that we actually plot (to avoid boundary effects)
# --- interstitial
mi = AtmModel(R = R, defect=Interstitial, ξ = ξi_off)
ui = solve(mi)
Dui = fdiff(ui, mi.B)
ri = bonddist(mi.X, mi.B)
I = find(ri .<= Rp)
ri = ri[I]; Dui = Dui[I]
# --- screw
ξ_s = ξs_off # [0.2; sqrt(3)/11]
ms = AtmModel(R = R, defect=Screw, ξ = ξ_s) # [0.3; sqrt(3)/12]
us = solve(ms)
Dus = fdiff(us, ms.B)
rs = bonddist(ms.X, ms.B)
I = find(rs .<= Rp)
rs = rs[I]; Dus = Dus[I]; es = ms.e[I]
# ---- plot
# to make this a nice plot, it takes a little bit of work
xbins = logspace(0, log10(Rp*0.99), 15)
xi, yi = MsUtils.envelope(ri, abs(Dui), xbins)
xs, ys = MsUtils.envelope(rs, abs(Dus), xbins)
xstot, ystot = MsUtils.envelope(rs, abs(es + Dus), xbins)
p = Axis([
Plots.Linear(xi, yi, legendentry="interstitial",
style="thick, blue, mark=*, mark options={fill=blue}");
Plots.Linear(xs, ys, legendentry="screw, corrector",
style="thick, red!70!black, mark=square*, mark options={fill=red!70!black}");
Plots.Linear(xstot, ystot, legendentry="screw, total",
style="thick, green!70!black, mark=triangle*, mark options = {fill=green!70!black}" );
MsUtils.plot_slope(6.0, 20.0, 0.5, -1.0);
Plots.Node(L"\sim r_b^{-1}", 12.0, 1e-1);
MsUtils.plot_slope(6.0, 20.0, 0.07, -2.0);
Plots.Node(L"\sim r_b^{-2}", 10.0, 2e-3);
],
style="xtick={1, 2, 4, 8, 16}, xticklabels={1,2,4,8,16}, ytickten={-1,-2,-3,-4}",
ymode="log", xmode="log",
xlabel=L"$r_b$",
ylabel = "strain (envelope)",
legendPos="south west" )
display(p); # save("ex-decay-all.pdf", p)
To avoid having to implement the Green's function, we implement the version of the LIN model described in Exercise 4.3.4, restricting displacements to a finite computational domain $$ \Omega_S = \Lambda \cap B_S $$ while the energy $E$ is approximated by $$ E^{\rm lin}(u) = \sum_{b \in \mathcal{B}^{\rm at}} \phi(e_b+Du_b) - \phi(e_b)
+ \sum_{b \in \mathcal{B} \setminus \mathcal{B}^{\rm at}}
\phi^{\rm lin}(e_b+Du_b) - \phi^{\rm lin}(e_b),
\phi^{\rm lin}(r) = \phi(0) + \phi'(0) r + \frac12 \phi''(0) r^2.
$$
In [ ]:
# Implementation of the LIN Model
# [1] Lin Model type : ------------------------------------------------
#
type LinModel{T}
R # domain parameter
S # domain parameter
X # lattice coordinates
Bat # nonlinear bonds
Blin # quadratic bonds
eat # pre-strains
elin # . . .
Ifree # free dofs
defect::T # type of defect
end
function LinModel(; R=5, S = 0.5 * R^3, defect = NoDefect, ξ=defcore(defect) )
# atoms
X = domain(S, defect; ξ=ξ)
# bonds
B = bonds(X)
rB = bonddist(X, B)
Iat = find(rB .<= R); Bat = (B[1][Iat], B[2][Iat])
Ilin = find(rB .> R); Blin = (B[1][Ilin], B[2][Ilin])
# prestrain
eat = prestrain(X, Bat, defect)
elin = prestrain(X, Blin, defect)
# return model
return LinModel(R, S, X, Bat, Blin, eat, elin, Ifree(X, S-1.1), defect)
end
# [2] energy
const ϕ1 = Dϕ(0)
const ϕ2 = D²ϕ(0)
ϕlin(r) = ϕ1 * r + 0.5 * ϕ2 * r.^2
Dϕlin(r) = ϕ1 + ϕ2 * r
D²ϕlin(r) = ϕ2 * ones(r)
energy(m::LinModel, U) = ( sum_kbn( ϕ(m.eat + fdiff(U[:], m.Bat)) - ϕ(m.eat) )
+ sum_kbn( ϕlin(m.elin + fdiff(U[:], m.Blin)) - ϕlin(m.elin) ) )
gradient(m::LinModel, U) = ( _grad_( m.Bat, Dϕ(m.eat + fdiff(U[:], m.Bat)), length(U) )
+ _grad_( m.Blin, Dϕlin(m.elin + fdiff(U[:], m.Blin)), length(U) ) )
hessian(m::LinModel, U) = ( _hess_( m.Bat, D²ϕ(m.eat + fdiff(U[:], m.Bat)), length(U) )
+ _hess_( m.Blin, D²ϕlin(m.elin + fdiff(U[:], m.Blin)), length(U) ) )
;
In [ ]:
# finite-difference check and Newton convergence check
m = LinModel(R = 4, S = 10, defect = Interstitial)
println("Finite-difference test of LIN Model")
fdtest(0.1 * rand(nsites(m)), x->energy(m,x), x->gradient(m,x), x->hessian(m,x))
println("\ntest Newton convergence")
ulin = solve(m, show=true);
Finally, we implement an atomistic/continuum coupling method. We give a very brief formulation of the model, but see the book manuscript, Sec. 4.5, for more details.
Displacements are now atomistic in a core region $\Lambda^{\rm at}$ but are continuous piecewise affine w.r.t. some triangulation $\mathcal{T}_h$ outside of the core region. The inner-most finite element nodes must lie on the lattice $\Lambda^{\rm at}$. Then we define atomistic bonds $\mathcal{B}^{\rm a}$ to be bonds between atoms in $\Lambda^{\rm at}$ without its boundary, interface bonds $\mathcal{B}^{\rm i}$ between interface atoms. The a/c energy is then given by $$ E^{\rm ac}h(u) = \sum{b \in \mathcal{B}^{\rm at}} \big(\phi(e_b + Du_b) - \phi(e_b)\big)
+ \sum_{b \in \mathcal{B}^{\rm i}} \frac12 \big(\phi(e_b + Du_b) - \phi(e_b)\big)
+ \sum_{T \in \mathcal{T}_h} |T| \big( W(E_T + \nabla u|_T) - W(E_T) \big),
$$ where $E_T$ denotes a continuum pre-strain field. In the interstitial case, $E_T = 0$ while in the screw dislocation case, $E_T = \nabla \hat{u}(x_T)$ where $\hat{u}$ is the predictor displacement.
This model now has three approximation parameters: the inner radius $R = R^{\rm at}$, the outer radius $R^{\rm c}$ and the mesh-size function $h(x)$. Both $R^{\rm c}$ and $h(x)$ are derived from $R^{\rm at}$ through error balancing.
In [ ]:
# Implementation of the A/C Scheme:
# our carefully constructed abstractions don't work so well directly;
# it is more convenient to re-implement some parts from scratch
# [1] pre-strain, continuum version
# recall the predicor: arg(x, y) = angle(x + im * y)
# the gradient is given by (-y,x) / |(x,y)|^2
"continuum pre-strain field for a screw dislocation"
prestrain_screw_cb(X) = [- X[2,:]; X[1,:]] ./ sumabs2(X,1) / (2π)
# [2] Model
"""
Atomistic/continuum coupling model
"""
type AcModel
R # atomistic region radius
Rc # continuum region radius
X # set of nodes (at + c)
tri # continuum region triangulation
Bat # atomistic bonds
Bi # interface bonds
eat # atomistic pre-strain
ei # interface pre-strain
Ec # continuum pre-strain
Ifree # free nodes
end
"""provide atomistic core and interface bonds separated"""
function acbonds(Tat)
# compute the atomistic bonds
Eat = [ Tat[1,:] Tat[1,:] Tat[2,:]; Tat[2,:] Tat[3,:] Tat[3,:] ]
Eat = sortcols(sort(Eat, 1))
occ = ones(Int, size(Eat, 2))
i = 1
while i < size(Eat, 2)
if Eat[:,i] != Eat[:,i+1]
occ[i] = 1
i += 1
else
occ[i] = 2
occ[i+1] = 0
i += 2
end
end
Bat = mat2tup(Eat[:, find(occ .== 2)])
Bi = mat2tup(Eat[:, find(occ .== 1)])
return Bat, Bi
end
function AcModel(; R=5, Rc=0.5*R^2, defect=NoDefect, ξ=defcore(defect))
Ra=R
# create the atomistic nodes
Xat = domain(Ra, defect; ξ=ξ)
# FEM nodes: reuse the implementation in LujiaLt
Xc, _ = LujiaLt.radial_nodes(Ra+0.4, Rc; hgrowth=1.0)
# concatenate
X = [Xat Xc]
# triangulate
tri = FEM.Triangulation(X)
# we could just get Bat from bonds(Xat), but then we
# will have difficulty deciding which bonds are core at
# and which are interface Instead, we use the triangulation
# compute all the information needed.
is_c = ones(FEM.nT(tri))
for (el, idx) in zip(FEM.elements(tri), 1:FEM.nT(tri))
x = tri.X[:, el.t]
# if all indices of the element are in Xat
# and if all edges are <= 1, then this is is an atomistic element
if (maximum(el.t) <= size(Xat,2)) &&
( (maximum(sumabs2(x[:,[2;3;1]] - x, 1)) <= 1.01) ||
(minimum(dist(x)) < 0.01) )
is_c[idx] = 0
end
end
# compute the atomistic bonds from this
Bat, Bi = acbonds(tri.T[:, find(is_c .== 0)])
# now get rid of the atomistic elements all-together
tri.T = tri.T[:, find(is_c .== 1)]
# precompute continuum strain
Ec = zeros(2, FEM.nT(tri))
if defect == Screw
for el in FEM.elements(tri)
Ec[:, el.idx] = prestrain_screw_cb(el.xT)
end
end
return AcModel( Ra, Rc, X, tri, Bat, Bi,
prestrain(X, Bat, defect), prestrain(X, Bi, defect), Ec,
Ifree(X, Rc-0.1) )
end
;
In [ ]:
include("MsUtils.jl")
# We can check-out what the A/C geometry looks like.
# this gives us some confidence that the implementation is correct
ac = AcModel(R=3.1, Rc=7, defect=Interstitial)
MsUtils.plot_acgeom(ac, plotwidth=10cm; ξ=ξi);
In [ ]:
# Implementation of the A/C Energy
# [3] First the pure Cauchy--Born part
# Cauchy Born potential
const Ann = ( [1.0;0.0], [0.5; sqrt(3)/2], [-0.5; sqrt(3)/2] )
const vv = 2 / sqrt(3)
W(F) = vv * sum( [ϕ(dot(F, a)) for a in Ann] )
DW(F) = vv * sum( [Dϕ(dot(F,a)) * a for a in Ann] )
D²W(F) = vv * sum( [D²ϕ(dot(F,a)) * a * a' for a in Ann] )
# compute the displacement gradient
∇u(el, U) = el.B' * U[el.t]
# we can use a comprehension for the energy . . .
energy_cb(m::AcModel, U) = sum_kbn(
Float64[ el.vol * (W(m.Ec[:, el.idx] + ∇u(el, U)) - W(m.Ec[:, el.idx]))
for el in FEM.elements(m.tri) ] )
# . . . but it is more difficult (not impossible!)
# to do the same for gradient and hessian
function gradient_cb(m::AcModel, U)
dE = zeros(U)
for el in FEM.elements(m.tri)
dE[el.t] += el.vol * el.B * DW(m.Ec[:, el.idx] + ∇u(el, U))
end
return dE
end
function hessian_cb(m::AcModel, U)
# initialise a triplet format
I = Int[]; J = Int[]; Z = Float64[];
for el in FEM.elements(m.tri)
# compute the local stiffness matrix
Aloc = el.vol * el.B * D²W(m.Ec[:, el.idx] + ∇u(el, U)) * el.B'
# and write it into the global matrix
for i = 1:3, j = 1:3
push!(I, el.t[i]); push!(J, el.t[j]); push!(Z, Aloc[i,j])
end
end
return sparse(I, J, Z, length(U), length(U))
end
# [4] total energy
energy(m::AcModel, U) =
sum_kbn( ϕ(m.eat + fdiff(U, m.Bat)) - ϕ(m.eat) ) +
0.5 * sum_kbn( ϕ(m.ei + fdiff(U, m.Bi)) - ϕ(m.ei) ) +
energy_cb(m, U)
gradient(m::AcModel, U) =
_grad_( m.Bat, Dϕ(m.eat + fdiff(U[:], m.Bat)), length(U) ) +
0.5 * _grad_( m.Bi, Dϕ(m.ei + fdiff(U[:], m.Bi)), length(U) ) +
gradient_cb(m, U)
hessian(m::AcModel, U) =
_hess_( m.Bat, D²ϕ(m.eat + fdiff(U[:], m.Bat)), length(U) ) +
0.5 * _hess_( m.Bi, D²ϕ(m.ei + fdiff(U[:], m.Bi)), length(U) ) +
hessian_cb(m, U)
;
In [ ]:
# finite-difference check and Newton convergence check
m = AcModel(R = 3, Rc = 6, defect = Interstitial)
println("Finite-difference test of AC Model")
fdtest(0.1 * rand(nsites(m)), x->energy(m,x), x->gradient(m,x), x->hessian(m,x))
println("\ntest Newton convergence")
uac = solve(m, show=true);
It remains to establish numerically that our predicted convergence rates are sharp. To this end, we implement some routines that allow us to estimate the distance (in energy-norm) between two computed solutions, $$ \big\| Du - Du_{\rm ex} \big\|_{\ell^2}, $$ possibly using different models and different computational domains. Here, $u_{\rm ex}$ is a high-accuracy solution to replace the unavailable exact solution.
Note that when estimating the error of a/c methods we will deviate slightly from the rigorous analysis where we have used $\| \cdot \|_{\rm ac}$ to measure the error.
In [ ]:
"""
`extend(Xlge, X, U)`
extend an approximate solution to be defined
on the nodes of a larger domain defined through m
"""
function extend(Xlge, X, U)
tri = FEM.Triangulation(X)
idx = FEM.locate(Xlge, tri)
Uext = zeros(length(idx))
for n = 1:length(idx)
if idx[n] == 0
Uext[n] = 0.0
else
t = tri.T[:, idx[n]]
lam = FEM.convex_coordinates(Xlge[:,n], idx[n], tri)
Uext[n] = dot(U[t], lam)
end
end
return Uext
end
"""
`compute_error(m, m_ex, U_ex, E_ex)`
compute the Ḣ¹ error and energy error of the model `m`,
as measured against a computation on a much larger domain
(m_ex, U_ex, E_ex).
"""
function compute_error(m, m_ex, U_ex, E_ex)
u = solve(m)
E = energy(m, u)
U_err = U_ex - extend(m_ex.X, m.X, u)
return norm( fdiff(U_err, m_ex.B) ), abs(E_ex - E)
end
Rout(m::AtmModel) = m.R
Rout(m::LinModel) = m.S
Rout(m::AcModel) = m.Rc
id(::AtmModel) = "ATM"
id(::LinModel) = "LIN"
id(::AcModel) = "AC"
"""
compute many errors for multiple models.
"""
function compute_errors(Model, Rlist, defect, ξ; show=false)
mlist = [ Model(R=R, defect=defect, ξ=ξ) for R in Rlist ]
# compute the exact solution
println("> Solving for 'exact' Solution ... ")
R_ex = 2 * maximum( [Rout(m) for m in mlist] )
print(" R_ex = $R_ex, "); sleep(0.1)
m_ex = AtmModel( R=R_ex, defect=defect, ξ=ξ )
print("#at = $(nsites(m_ex)) ... ")
show && println()
u_ex = solve(m_ex, show=show)
nrm_ex = norm(fdiff(u_ex, m_ex.B))
E_ex = energy(m_ex, u_ex)
println("DONE.")
# now loop through approximations
err2 = Float64[]
errE = Float64[]
for m = mlist
print("> Solving for $(id(m)), #at = $(nsites(m)) ... "); sleep(0.1)
err2_, errE_ = compute_error(m, m_ex, u_ex, E_ex)
push!(err2, err2_/nrm_ex)
push!(errE, errE_/abs(E_ex))
println("DONE.")
end
return err2, errE
end
;
In [ ]:
#### INTERSTITIAL TEST
defect = Interstitial
ξ = [0.4; 0.1]
# atm test
R_atm = [4; 6; 8; 10; 13; 16; 20; 25]+0.1
err2_atm, errE_atm = compute_errors(AtmModel, R_atm, defect, ξ)
# lin test
R_lin = [2; 3; 4; 5]+0.1
err2_lin, errE_lin = compute_errors(LinModel, R_lin, defect, ξ)
# ac test
R_ac = [2; 3; 4; 6; 8; 10]+0.1
err2_ac, errE_ac = compute_errors(AcModel, R_ac, defect, ξ)
;
In [ ]:
####### PLOT OF ENERGY-NORM ERROR / INTERSTITIAL
using PGFPlots
myblue="blue"; myred="red!70!black"; mygreen="green!70!black";
p = Axis([
Plots.Linear(R_atm, err2_atm, legendentry="ATM",
style="thick, $myblue, mark=*, mark options={draw=black,fill=$myblue}");
Plots.Linear(R_lin, err2_lin, legendentry="LIN",
style="thick, $myred, mark=square*, mark options={fill=$myred}");
Plots.Linear(R_ac, err2_ac, legendentry="AC",
style="thick, $mygreen, mark=triangle*, mark options = {fill=$mygreen}" );
MsUtils.plot_slope(10.0, 22.0, 0.22, -1.0);
Plots.Node(L"\sim R^{-1}", 15.0, 9e-3);
MsUtils.plot_slope(5.5, 9.0, 0.4, -2.0);
Plots.Node(L"\sim R^{-2}", 7.0, 4e-3);
MsUtils.plot_slope(2.8, 4.2, 0.3, -3.0);
Plots.Node(L"\sim R^{-3}", 2.5, 6e-3);
],
style="""xtick={2,3,5,10,17,25}, xticklabels={2,3,5,10,17,25},
ytickten={-2.5, -2, -1.5, -1}, minor ytick={0.001}""",
ymode="log", xmode="log",
xlabel=L"$R$",
ylabel = L"\|D\bar{u} - D\bar{u}^{\rm apx} \|_{\ell^2}",
title="Interstitial : Energy-Norm Errors",
legendPos="north east" )
display(p); #save("apx_err2_int.pdf", p)
In [ ]:
####### PLOT OF ENERGY ERROR / INTERSTITIAL
p = Axis([
Plots.Linear(R_atm, errE_atm, legendentry="ATM",
style="thick, $myblue, mark=*, mark options={draw=black,fill=$myblue}");
Plots.Linear(R_lin, errE_lin, legendentry="LIN",
style="thick, $myred, mark=square*, mark options={fill=$myred}");
Plots.Linear(R_ac, errE_ac, legendentry="AC",
style="thick, $mygreen, mark=triangle*, mark options = {fill=$mygreen}" );
MsUtils.plot_slope(10.0, 22.0, 0.22, -2.0);
Plots.Node(L"\sim R^{-2}", 18.0, 1.5e-3);
MsUtils.plot_slope(6.0, 11.0, 0.9, -4);
Plots.Node(L"\sim R^{-4}", 12, 2e-4);
MsUtils.plot_slope(2.9, 4.4, 0.18, -6.0);
Plots.Node(L"\sim R^{-6}", 2.8, 4e-5);
],
style="""xtick={2,3,5,10,17,25}, xticklabels={2,3,5,10,17,25},
ytickten={-4.8,-4,-3,-2}, minor ytick={0.001}""",
ymode="log", xmode="log",
xlabel=L"$R$",
ylabel = L"|E(\bar{u}) - E^{\rm apx}(\bar{u}^{\rm apx}) |",
title="Interstitial : Energy Errors",
legendPos="north east" )
display(p); # save("apx_errE_int.pdf", p)
In [ ]:
#### SCREW TEST
defect = Screw
ξ = [0.25; 0.12]
# atm test
R_atm = [8; 10; 13; 16; 20; 25]+0.1
err2_atm, errE_atm = compute_errors(AtmModel, R_atm, defect, ξ)
# lin test
R_lin = [3; 4; 5; 6; 7]+0.1
err2_lin, errE_lin = compute_errors(LinModel, R_lin, defect, ξ.+0.1)
# ac test
R_ac = [4; 6; 8; 10; 13; 16]+0.1
err2_ac, errE_ac = compute_errors(AcModel, R_ac, defect, ξ)
;
In [ ]:
####### PLOT OF ENERGY-NORM ERROR / SCREW
myblue="blue"; myred="red!70!black"; mygreen="green!70!black";
p = Axis([
Plots.Linear(R_atm, err2_atm, legendentry="ATM",
style="thick, $myblue, mark=*, mark options={draw=black,fill=$myblue}");
Plots.Linear(R_lin, err2_lin, legendentry="LIN",
style="thick, $myred, mark=square*, mark options={fill=$myred}");
Plots.Linear(R_ac, err2_ac, legendentry="AC",
style="thick, $mygreen, mark=triangle*, mark options = {fill=$mygreen}" );
MsUtils.plot_slope(10.0, 22.0, 1.0, -1.0);
Plots.Node(L"\sim R^{-1}", 18.0, 8.0e-2);
MsUtils.plot_slope(9.0, 17.0, 2.4, -2.0);
Plots.Node(L"\sim R^{-2}", 15.5, 1.8e-2);
MsUtils.plot_slope(2.8, 3.5, 1.0, -3.0);
Plots.Node(L"\sim R^{-3}", 3.0, 2e-2);
],
style="""xtick={2,3,5,10,17,25}, xticklabels={2,3,5,10,17,25},
ytickten={-2.5, -2, -1.5, -1}, minor ytick={0.001}""",
ymode="log", xmode="log",
xlabel=L"$R$",
ylabel = L"\|D\bar{u} - D\bar{u}^{\rm apx} \|_{\ell^2}",
title="Screw Dislocation : Energy-Norm Errors",
legendPos="south west" )
display(p); # save("apx_err2_screw.pdf", p)
In [ ]:
####### PLOT OF ENERGY ERROR / SCREW
p = Axis([
Plots.Linear(R_atm, errE_atm, legendentry="ATM",
style="thick, $myblue, mark=*, mark options={draw=black,fill=$myblue}");
Plots.Linear(R_lin, errE_lin, legendentry="LIN",
style="thick, $myred, mark=square*, mark options={fill=$myred}");
Plots.Linear(R_ac[1:end-1], errE_ac[1:end-1], legendentry="AC",
style="thick, $mygreen, mark=triangle*, mark options = {fill=$mygreen}");
MsUtils.plot_slope(12.0, 25.0, 0.8, -2.0);
Plots.Node(L"\sim R^{-2}", 17.0, 5.8e-3);
MsUtils.plot_slope(6.3, 11.5, 0.6, -4);
Plots.Node(L"\sim R^{-4}", 6, 1.6e-4);
#MsUtils.plot_slope(2.9, 4.4, 0.18, -6.0);
#Plots.Node(L"\sim R^{-3}", 2.8, 4e-5);
],
style="""xtick={2,3,5,10,17,25}, xticklabels={2,3,5,10,17,25},
ytickten={-4.8,-4,-3,-2}, minor ytick={0.001}""",
ymode="log", xmode="log",
xlabel=L"$R$",
ylabel = L"|E(\bar{u}) - E^{\rm apx}(\bar{u}^{\rm apx}) |",
title="Screw Dislocation : Energy Errors",
legendPos="south east" )
display(p); # save("apx_errE_screw.pdf", p)