From 88daf3e1d58a1041b7c24905fe2b65768eed0fb0 Mon Sep 17 00:00:00 2001 From: Stephan Hilb <stephan@ecshi.net> Date: Sat, 26 Feb 2022 13:08:20 +0100 Subject: [PATCH] fixup error estimator slightly --- scripts/run_experiments.jl | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl index 62ecfe4..03f35a1 100644 --- a/scripts/run_experiments.jl +++ b/scripts/run_experiments.jl @@ -622,7 +622,6 @@ function estimate_res!(st::L1L2TVState) st.T(adjoint, tdata, p1) + (isSnabla ? zero(p1) : st.beta * u) facetf(hfacet, n; g, u, nablau, p1, p2, tdata) = - #display(n' * p2) # norm2 is calculated later after both facet contributions are # consolidated n' * ((isSnabla ? st.beta * st.S(zero(u), nablau) : zero(p2)) + p2) @@ -633,7 +632,8 @@ function estimate_res!(st::L1L2TVState) fs = (;st.g, st.u, nablau = nabla(st.u), st.p1, st.p2, st.tdata) space = st.est.space - facetV = Dict{Tuple{Int, Int}, MVector{ndims_u_codomain(st), Float64}}() + facetv = Dict{Tuple{Int, Int}, MVector{ndims_u_codomain(st), Float64}}() + st.est.data .= 0. for cell in cells(mesh) A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}( @@ -670,29 +670,33 @@ function estimate_res!(st::L1L2TVState) v = A[:, j] - A[:, i] hfacet = norm(v) normal = -normalize(cross(v)) - intel = hfacet - # we don't need to recompute values since only piecewise + # we don't need to recompute `opvalues` since only piecewise # constant data is used in `facetf` - # we use sqrt here since this value will be squared later - facetres = (isSnabla ? sqrt(hfacet) : inv(sqrt(hfacet))) * - facetf(hfacet, normal; opvalues...) .* sqrt(intel) + facetres = facetf(hfacet, normal; opvalues...) # average facet contributions - v = SVector(facetres) - facetV[(fi, fj)] = get(facetV, (fi, fj), zero(v)) + SVector(v) - facetV[(fj, fi)] = get(facetV, (fj, fi), zero(v)) + SVector(v) + v = SVector{ndims_u_codomain(st)}(facetres / 2) + facetv[(fi, fj)] = get(facetv, (fi, fj), zero(v)) + v + facetv[(fj, fi)] = get(facetv, (fj, fi), zero(v)) + v end end # add facet contributions to cells for cell in cells(mesh) for (i, j) in ((i, mod1(i + 1, 3)) for i in 1:3) + A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}( + view(mesh.vertices, :, view(mesh.cells, :, cell))) + v = A[:, j] - A[:, i] + hfacet = norm(v) + intel = hfacet + fi = mesh.cells[i, cell] fj = mesh.cells[j, cell] gdofs = space.dofmap[:, 1, cell] - st.est.data[gdofs] .+= norm2(facetV[(fi, fj)]) + st.est.data[gdofs] .+= (isSnabla ? hfacet : inv(hfacet)) * + norm2(facetv[(fi, fj)]) * intel end end -- GitLab