Select Git revision
check-package.m4
run.jl 3.77 KiB
export myrun
using LinearAlgebra: norm
function myrun()
name = "test"
mesh = init_grid(50)
d = 2
m = 1
Vg = FeSpace(mesh, DP0(), (1,))
Vu = FeSpace(mesh, P1(), (m,))
Vp1 = FeSpace(mesh, DP0(), (1,))
Vp2 = FeSpace(mesh, DP1(), (m, d))
g = FeFunction(Vg, name="g")
mask = FeFunction(Vg, name="mask")
u = FeFunction(Vu, name="u")
p1 = FeFunction(Vp1, name="p1")
p2 = FeFunction(Vp2, name="p2")
du = FeFunction(Vu)
dp1 = FeFunction(Vp1)
dp2 = FeFunction(Vp2)
nablau = nabla(u)
nabladu = nabla(du)
g.data .= 0
u.data .= 0
p1.data .= 0
p2.data .= 0
du.data .= 0
dp1.data .= 0
dp2.data .= 0
alpha1 = 0.
alpha2 = 10.
beta = 0.
lambda = 0.01
gamma1 = 1e-3
gamma2 = 1e-3
interpolate!(g, x -> norm(x - SA[0.5, 0.5]) < 0.3)
interpolate!(mask, x -> abs(x[2] - 0.5) > 0.1)
#T(data, u) = u
#tdata = u # FIXME
@inline T(data, u) = iszero(data) ? zero(u) : u
tdata = mask
#T(data, u) = dot(data, u)
#tdata = nabla(fw)
S(u) = u
function dp1_update(x; g, u, p1, du, tdata)
m1 = max(gamma1, norm(T(tdata, u) - g))
cond = norm(T(tdata, u) - g) > gamma1 ?
dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
zeros(size(p1))
return -p1 + alpha1 / m1 * (T(tdata, u) + T(tdata, du) - g) - cond
end
function dp2_update(x; u, nablau, p2, du, nabladu)
m2 = max(gamma2, norm(nablau))
cond = norm(nablau) > gamma2 ?
dot(nablau, nabladu) / norm(nablau)^2 * p2 :
zeros(size(p2))
return -p2 + lambda / m2 * (nablau + nabladu) - cond
end
function p1_project!()
p1.space.element::DP0
p1.data .= clamp.(p1.data, -alpha1, alpha1)
end
function p2_project!()
p2.space.element::DP1
p2d = reshape(p2.data, m * d * 3, :) # no copy
for i in axes(p2d, 2)
p2in = norm(p2d[:, i])
if p2in > lambda
p2d[:, i] .*= lambda ./ p2in
end
end
end
@inline function du_a(x, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2, tdata)
m1 = max(gamma1, norm(T(tdata, u) - g))
cond1 = norm(T(tdata, u) - g) > gamma1 ?
dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
zeros(size(p1))
a1 = alpha1 / m1 * dot(T(tdata, du), T(tdata, phi)) -
dot(cond1, T(tdata, phi))
m2 = max(gamma2, norm(nablau))
cond2 = norm(nablau) > gamma2 ?
dot(nablau, nabladu) / norm(nablau)^2 * p2 :
zeros(size(p2))
a2 = lambda / m2 * dot(nabladu, nablaphi) -
dot(cond2, nablaphi)
aB = alpha2 * dot(T(tdata, du), T(tdata, phi)) +
beta * dot(S(du), S(phi))
return a1 + a2 + aB
end
@inline function du_l(x, phi, nablaphi; g, u, nablau, tdata)
aB = alpha2 * dot(T(tdata, u), T(tdata, phi)) +
beta * dot(S(u), S(phi))
m1 = max(gamma1, norm(T(tdata, u) - g))
p1part = alpha1 / m1 * dot(T(tdata, u) - g, T(tdata, phi))
m2 = max(gamma2, norm(nablau))
p2part = lambda / m2 * dot(nablau, nablaphi)
gpart = alpha2 * dot(g, T(tdata, phi))
return -aB - p1part - p2part + gpart
end
function save_step!(pvd, i)
vtk = vtk_mesh("$(name)_$(lpad(i, 5, '0')).vtu", mesh)
vtk_append!(vtk, g, u, p1, p2)
vtk_save(vtk)
pvd[i] = vtk
end
pvd = paraview_collection("$(name).pvd")
save_step!(pvd, 0)
for i = 1:20
print("[$(lpad(i, 3, '0'))] ")
# solve du
print("assemble ... ")
A = assemble(Vu, du_a; g, u, nablau, p1, p2, tdata)
b = assemble_rhs(Vu, du_l; g, u, nablau, tdata)
print("solve ... ")
du.data .= A \ b
# solve dp1, dp2
interpolate!(dp1, dp1_update; g, u, p1, du, tdata)
interpolate!(dp2, dp2_update; u, nablau, p2, du, nabladu)
# newton update
u.data .+= du.data
p1.data .+= dp1.data
p2.data .+= dp2.data
# reproject p1, p2
p1_project!()
p2_project!()
print("save ... ")
save_step!(pvd, i)
println()
end
vtk_save(pvd)
return
end