From 1e93e03adff07cdb421381f078d397d6468f870c Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Mon, 28 Feb 2022 19:43:55 +0100
Subject: [PATCH] fix nablau dimension issue

---
 scripts/run_experiments.jl | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 9b116f4..8a316a2 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -459,7 +459,7 @@ function step_pd2!(st::L1L2TVState; sigma, tau, theta = 1.)
 
     function p2_update(x_; p2, nablau)
 	st.lambda == 0 && return zero(p2)
-	return (p2 + sigma * nablau) /
+	return (p2 + sigma * nablau') /
 	    (1 + st.gamma2 * sigma / st.lambda)
     end
     interpolate!(p2_next, p2_update; st.p2, nablau = nabla(u_new))
@@ -500,7 +500,7 @@ function step_pd2!(st::L1L2TVState; sigma, tau, theta = 1.)
     u_lp(x, phi, nablaphi; g, u, p1, p2, tdata) =
 	- tau * (
 	    dot(p1, st.T(tdata, phi)) +
-	    dot(p2, nablaphi))
+	    dot(p2, nablaphi'))
 
     Au, bu = assemble(u_new.space, u_au, u_lu; st.g, st.u, st.p1, st.p2, st.tdata)
     Ap, bp = assemble_lumped(u_new.space, u_ap, u_lp; st.g, st.u, st.p1, st.p2, st.tdata)
@@ -610,17 +610,16 @@ end
 
 # this computes the residual error indicator which has missing theoretical
 # support for alpha1 > 0
-# TODO: finish!
 function estimate_res!(st::L1L2TVState)
     # FIXME: dirty hack
-    isSnabla = st.S(0, 1) == 1
+    isSnabla = st.S(0, 1) == 1'
     norm2(v) = dot(v, v)
 
     cellf(hcell; g, u, nablau, p1, p2, tdata) =
          st.alpha2 *
             st.T(adjoint, tdata, st.alpha2 * st.T(tdata, u) - g) +
             st.T(adjoint, tdata, p1) +
-            (isSnabla ? zero(p1) : st.beta * u)
+            (isSnabla ? zero(u) : st.beta * u)
     facetf(hfacet, n; g, u, nablau, p1, p2, tdata) =
         # norm2 is calculated later after both facet contributions are
         # consolidated
@@ -661,7 +660,6 @@ function estimate_res!(st::L1L2TVState)
 	st.est.data[gdofs] .= cellres
 
         # facet contributions
-
         cross(v) = SA[-v[2], v[1]]
         for (i, j) in ((i, mod1(i + 1, 3)) for i in 1:3)
             fi = mesh.cells[i, cell]
-- 
GitLab