

In [ ]:

    
using LujiaLt
using LujiaLt.Potentials, LujiaLt.Plotting, LujiaLt.Solve
using PyPlot, Plots, Compose

"""
generate Atm model with domain with radius R. The `edgevacancy` parameter 
stabilises the edge dislocation so that the computations are easier and more robust
"""
EdgeDislocation(R; kwargs...) = 
    Atm(V=LennardJonesPotential(;kwargs...), Ra = R, 
        defect=:edge, edgevacancy=true)
;

Test geometry construction



In [ ]:

    
# Test setup of an edge dislocation configuration via LujiaLt
at = EdgeDislocation(8.1)
LujiaLt.Plotting.plot(at, plotwidth=12cm, bondwidth=0.3)

Testing force decay and net-force decay

Checking that $|F(\ell)| \lesssim |\ell|^{-3}$ and $\sum_{\ell \in \Lambda} F(\ell) = 0$, which are the requirements for the convergence theory.



In [ ]:

    
"""
`forcedecay(R::Real) -> (r, f)`

where `r` is a vector of distances of atoms from the core and `f` 
a vector containing the magnitude of the forces acting on the atoms. 
`r` and `f` are only returned for atoms *inside* the domain, 
boundary atoms are excluded.
"""
function forcedecay(R::Real; kwargs...)
    at = EdgeDislocation(R; kwargs...)
    F = forces(at)
    f = sumabs2(F, 1) |> sqrt
    r = sumabs2(at.Yref, 1) |> sqrt
    I = find(r .<= R)
    return r[I], f[I]
end

"""
`netforce(R::Real; varargs...) -> Fnet`

`Fnet` is the sum over all forces, excluding the boundary forces.
"""
function netforce(R::Real; kwargs...)
    at = EdgeDislocation(R; kwargs...)
    return sum(forces(at), 2)
end 
;



In [ ]:

    
cutoff = (1.5, 2.5)    # (2.2, 2.8)

subplot(2,2,1)
RR = [2^p + 0.1 for p in 2:7 ]   # 8
nF = [ norm( netforce(R, cutoff=cutoff) ) for R in RR ]
PyPlot.loglog(RR, nF, "bo-", RR, 5*RR.^(-2), "r:")
ylabel("|net-force|")
xlabel("R - domain radius")
title("Convergence of net-Force")

subplot(2,2,2)
R = 99.1     # 399.1
r, f = forcedecay(R, cutoff=cutoff)
loglog(r, f, "b.", [1, R], 200*[1,R].^(-3), "r:", markersize=1)
legend((L"|F|", L"\sim r^{-3}"), loc="lower left")
xlabel("r")
ylabel("|F|")
axis([0.8, R*1.2, 1e-6, 1e2])
title("Decay of Forces")

Solver test

Quick check that the nonlinear solver works as expected.



In [ ]:

    
at = EdgeDislocation(10.1)
Y = solve(at, show_trace=false, display_result = true);
# If this fails, try to `Pkg.pin("Optim", v"0.6.1")`



In [ ]:

    
# note that the geometry is O(1) different from the CLE configuration
LujiaLt.Plotting.plot(at, X=Y, plotwidth=12cm, bondwidth=0.3)

Convergence Test

test the predicted convergence rate - this does not work yet. need to find the bug; first, check the correct decay of the corrector solution?



In [ ]:

    
# domain radii
RR = [2.1, 4.1, 8.1, 16.1, 32.1, 64.1]
Rex = 3 * maximum(RR)

# "exact" solution
print("Compute comparison solution ... "); sleep(0.01)
atex = EdgeDislocation(Rex)
Yex = solve(atex)
println("done.")

# loop through smaller systems 
err2 = zeros(length(RR))
for (i, R) in enumerate(RR)
    print("compute R = $R ... "); sleep(0.01)
    at = EdgeDislocation(R)
    Y = solve(at)
    err2[i] = LujiaLt.error_energynorm(Y, at, Yex, atex)
    println("done.")
end

# save computation
save("edgelj_errors.jld", "RR", RR, "err2", err2)



In [ ]:

    
using JLD
RR, err2 = load("edgelj_errors.jld", "RR", "err2")
dofs = 3 * (RR-0.1).^2 * 2
using Plots
Plots.pyplot()
fnt = Plots.font(14, "Times")
Plots.plot(dofs, err2, lw=2.0, marker = (:circle, 10), label = L"\|\nabla \bar{u} - \nabla \bar{u}_R\|_{L^2}", 
    xaxis = (L"$\#$ dofs", (dofs[1]/1.3, dofs[end]*1.3), :log, fnt), 
    yaxis = ("error", :log, (8e-3, 3.0), fnt),
    xticks = (dofs, [string(round(Int, d)) for d in dofs]), xtickfont = fnt, yticks = [-2, -1, 0], ytickfont = fnt )
Plots.plot!(dofs[3:5], 2*dofs[3:5].^(-1/2), line = (:dash, 1.0), label=L"$\sim R^{-1/2}$")
Plots.title!("Galerkin Approximation of Edge Dislocation", titlefont = fnt, legendfont = fnt)
Plots.savefig("edgelj_errors.pdf")



In [ ]: