EquilibratedFlux.jl
This package is based on Gridap.jl to provide post-processing tools to calculate reconstructed fluxes associated to the given approximate solution of a PDE.
For simplicity, we consider here the Poisson equation
\[\begin{align} - \Delta u &= f &&\text{in }\Omega\\ u &= g &&\text{on }\partial\Omega. \end{align}\]
We suppose we have already computed a conforming approximation $u_h \in V_h\subset H^1_0(\Omega)$ to the solution $u$ in Gridap.jl by solving
\[(\nabla u_h, \nabla v_h) = (f, v_h)\quad\forall v_h\in V_h,\]
The EquilibratedFlux.jl library then provides the tools to compute a reconstructed flux associated to $u_h$. This flux, obtained by postprocessing, is an approximation to the numerical flux, i.e.
\[\sigma_h \approx -\nabla u_h.\]
This flux has the important property of being "conservative over faces" in the sense that
\[\sigma_h \in \mathbf{H}(\mathrm{div},\Omega).\]
We provide two functions to obtain such an object: build_equilibrated_flux and build_averaged_flux both provide reconstructed fluxes, which we denote by $\sigma_{\mathrm{eq},h}$ and $\sigma_{\mathrm{ave},h}$ respectively.
In addition to the properties listed above, the equilibrated flux $\sigma_{\mathrm{eq},h}$ satisfies the so-called equilibrium condition, i.e., for piecewise polynomial $f$, we have
\[\nabla\cdot\sigma_{\mathrm{eq},h} = f.\]
The reconstructed flux is the main ingredient in computing a posteriori error estimators. See the first tutorial for a complete demonstration of how to do this.
EquilibratedFlux.build_averaged_flux — Method
build_averaged_flux(𝐀ₕ, model)Builds a simple lowest-order (RT₀) H(div, Ω)-conforming flux reconstruction σ_ave by averaging the (generally discontinuous, face-wise two-valued) normal component of 𝐀ₕ across each interior face, and taking it directly on boundary faces. 𝐀ₕ is typically -∇(uh) for a conforming approximation uh of the Poisson problem. Unlike build_equilibrated_flux, σ_ave does not satisfy the equilibrium property ∇⋅σ_ave = f, so the resulting a posteriori estimator ‖σ_ave - 𝐀ₕ‖ is typically less sharp; see the package tutorials for a comparison.
EquilibratedFlux.build_equilibrated_flux — Method
build_equilibrated_flux(𝐀ₕ, f, model, RT_order; measure = nothing, weight= 1.0,
neumann_tags = String[], neumann_data = nothing)Builds the equilibrated flux reconstruction σ_eq ∈ H(div, Ω) (as a Raviart-Thomas FEFunction of order RT_order) associated with a conforming approximation uh of the Poisson problem -∇⋅(weight*∇u) = f, by solving small mixed finite element problems on vertex-centered patches of model. 𝐀ₕ is the numerical flux to reconstruct, typically -∇(uh) (or -weight*∇(uh) for a non-unit weight); f is the source term. σ_eq satisfies ∇⋅σ_eq = f exactly for piecewise polynomial f of degree RT_order, and σ_eq ≈ 𝐀ₕ in the sense that ‖σ_eq - 𝐀ₕ‖ provides a reliable and efficient a posteriori error estimator for ‖∇(u - uh)‖ (see the package tutorials for details). measure overrides the default quadrature (degree = 2*RT_order + 2) used to assemble the patch problems.
neumann_tags marks the part ΓN of the boundary (as model face labeling tags) where an inhomogeneous Neumann condition -∇u⋅n = σN holds. When given, neumann_data must be a function x -> VectorValue giving the physical flux σN(x)⋅n(x) (e.g. -∇u_exact for a manufactured solution), matching the 𝐀ₕ = -∇(uh) convention used elsewhere in this package.
TODO: mixed Dirichlet/Neumann boundaries are not yet supported: neumann_tags must currently mark either the entire boundary (pure Neumann) or none of it (pure Dirichlet), and an error is raised otherwise. See the corresponding TODOs in Patch.jl (create_patches) for what interface-vertex support would require.