diff --git a/FD/src/FD.jl b/FD/src/FD.jl
index d5621fcb1b556f27d858ffe64d537fb79a011e25..794a1ad1867a253527992dd7d34e47a88cfade05 100644
--- a/FD/src/FD.jl
+++ b/FD/src/FD.jl
@@ -225,11 +225,19 @@ grid(f::GridFunction) = f.grid
Base.view(f::AbstractGridFunction{G}, idx...) where G<:MetaGrid =
GridFunction(GridView(f.grid, indices(f.grid, idx...)), view(f.data, indices(f.grid, idx...)..., :))
-Base.parent(f::AbstractGridFunction{<:GridView}) =
+function zero_continuation(f::GridFunction{<:GridView})
+ dst = GridFunction(parent(grid(f)), zeros(size(parent(grid(f)))..., rangedim(f)), Val(rangedim(f)))
+ dst.data[grid(f).indices..., :] .= f.data
+ dst
+end
+
+Base.parent(f::GridFunction{<:GridView,<:Any,<:Any,<:Any,<:SubArray}) =
GridFunction(parent(grid(f)), parent(f.data))
+
+
function _check_shape(a::GridFunction, b::GridFunction)
#TODO: why does fv.grid == gv.grid fail?
a.grid === b.grid || throw(ArgumentError("grids need to match"))
diff --git a/chambollefp_dd.jl b/chambollefp_dd.jl
index 10e73b0ce31eafe1d05747f2e77122e7ea1a827d..b622ff9c630d5140a8307159ebf32424b57c67aa 100644
--- a/chambollefp_dd.jl
+++ b/chambollefp_dd.jl
@@ -36,35 +36,14 @@ end
# TODO: investigate matrix multiplication performance with GridFunctions
function fpiter(grid, p, wg, A, α, β, τ)
- ## local variant
- #w = FD.GridFunction(grid, view(wg.data, grid.indices..., :))
- #x = FD.GridFunction(grid, (A'*A + β*I) \ (FD.divergence(p) + w).data[:])
- #q = FD.gradient(x)
-
- ## global variant
- x = FD.GridFunction(FD.parent(grid), (A'*A + β*I) \ (FD.divergence(p) + wg).data[:])
- #FD.my_plot(FD.divergence(p) + wg)
-
+ x = FD.GridFunction(FD.parent(grid), (A'*A + β*I) \ (FD.divergence_z(FD.zero_continuation(p)) + wg).data[:])
q = FD.GridFunction(grid, view(FD.gradient(x).data, grid.indices..., :))
- q2 = FD.gradient(FD.GridFunction(grid, view(x.data, grid.indices..., :)))
-
- #display(x.data[10:13,51:54,1])
- #display(q.data[10:13,52:53,2])
- #display(q2.data[10:13,52:53,2])
- #sleep(0.1)
-
+ ## this misses out on the right boundary
#q = FD.gradient(FD.GridFunction(grid, view(x.data, grid.indices..., :)))
-
- #FD.my_plot(FD.pw_norm(q))
- #FD.my_plot(x)
- #FD.my_plot(α)
-
p_new = FD.GridFunction(grid, FD.pwmul(FD.GridFunction(grid, α ./ (α .+ τ .* FD.pw_norm(q))),
FD.GridFunction(grid, p .+ τ .* q)))
-
-
reserr = maximum(sum((p_new .- p).^2, dims=ndims(p)))
(p_new, reserr)
end
@@ -72,113 +51,47 @@ end
function fpalg(grid, w, A, p = FD.GridFunction(x -> 0, ndims(grid), grid);
α, β, τ = 1/8, ϵ = 1e-4)
- local uimg
reserr = 1
k = 0
while reserr > ϵ && k < 10e1
k += 1
(p, reserr) = fpiter(grid, p, w, A, α, β, τ)
-
- #u = p2u(grid, p, A, f, α=α, β=β)
- #@show k
- #FD.my_plot(u)
end
return p
end
function dditer(grid, p, A, f; α, β, τ, ϵ, ρ)
- #pzero = FD.GridFunction(x -> 0, ndims(grid), grid)
-
p̂ = (1-ρ) * p
i = 0
for part in FD.parts(grid)
gridv = view(grid, Tuple(part)...)
+
fv = view(f, Tuple(part)...)
- #p.data .= 0
+
p̂v = view(p̂, Tuple(part)...)
pv = view(p, Tuple(part)...)
- #pzerov = view(pzero, Tuple(part)...)
θ = theta(grid, Tuple(part)...)
θv = view(θ, Tuple(part)...)
- #y = FD.divergence(FD.pwmul(FD.GridFunction(grid, 1 .- θ), p))
- #FD.my_plot(y)
- #sleep(1)
-
w = FD.GridFunction(grid, A' * f.data[:]) + FD.divergence(FD.pwmul(FD.GridFunction(grid, 1 .- θ), p))
wv = view(w, Tuple(part)...)
- #FD.my_plot(w)
- #sleep(1)
-
- #wv = A' * fv
-
- #u = p2u(grid, p̂, A, f, α=α, β=β)
- #uv = view(u, Tuple(part)...)
-
- #FD.my_plot(uv)
-
- #sleep(1)
-
p̂v .+= ρ .* fpalg(gridv, w, A, α=α*θv, β=β, τ=τ, ϵ=ϵ)
- #p̂v .= fpalg(gridv, wv, A, α=α*θv, β=β, τ=τ, ϵ=ϵ)
-
- #FD.my_plot(FD.pw_norm(p))
- #sleep(1)
- #FD.my_plot(FD.pw_norm(f))
- #sleep(1)
- #FD.my_plot(FD.pw_norm(p̂))
- #display(p̂.data[10:13,52:54,:])
- #break
- #FD.my_plot(FD.pw_norm(p̂))
- #sleep(1)
- #FD.my_plot(θ)
- #display(p̂v.data)
- #break
- #sleep(1)
-
- #u = p2u(grid, p, A, f, α=α*θ, β=β)
- u = p2u(grid, p̂, A, f, α=α, β=β)
- uv = view(u, Tuple(part)...)
-
- #FD.my_plot(θv)
-
- #sleep(1)
- #FD.my_plot(FD.pw_norm(p̂v))
- #@show maximum(u), minimum(u)
- #@show maximum(p̂v), minimum(p̂v)
- #sleep(1)
- #FD.my_plot(FD.pw_norm(p̂))
- #FD.my_plot(u)
- #FD.my_plot(p̂v)
- #y = FD.GridFunction(gridv, FD.gradient(θv).data[:,:,1])
- #display(y.data)
- #display(θv.data)
- #FD.my_plot(FD.GridFunction(gridv, FD.gradient(θv).data[:,:,2]))
-
- #display(FD.gradient(θv).data)
- #break
-
- #@show size(θv.data), size(p̂v.data)
-
- #FD.my_plot(FD.pw_norm(p̂))
i += 1
- #break
- #i > 3 && break
end
p .= p̂
u = p2u(grid, p, A, f, α=α, β=β)
- #FD.my_plot(u)
- FD.my_plot(FD.pw_norm(p))
+ FD.my_plot(u)
+ #FD.my_plot(FD.pw_norm(p))
end
-grid = FD.MetaGrid((0:0.01:1, 0:0.01:1), (1, 2), 5)
+grid = FD.MetaGrid((0:0.01:1, 0:0.01:1), (2, 2), 5)
A = I
#A = spdiagm(0 => ones(length(grid)), 1 => ones(length(grid)-1) )
@@ -195,7 +108,7 @@ g = FD.GridFunction(x -> norm(x), 2, grid)
p = FD.GridFunction(x -> 0, 2, grid)
-for n = 1:1
+for n = 1:10
println("iteration $n")
dditer(grid, p, A, f, α=α, β=β, τ=τ, ϵ=ϵ, ρ=ρ)
end