Introduce Flux Error Estimation

[hughcars]

Introduce flux error estimation based on a global ZZ solve.

• Introduces the FluxProjector class, which is used for computing a "coarse flux" in
  a vector L2 space, then performing a global mass matrix inversion to compute a
  "smooth" flux that better represents the underlying physical quantity. This is achieved
  through a multigrid hierarchy.
• Introduces the CurlFluxErrorEstimator and GradFluxErrorEstimator classes, which utilize the
  FluxProjector class to compute elemental error estimates comparing the coarse numerical flux
  from the solution against a smoothed flux from a global mass matrix inversion.
• Introduces the CurlFluxErrorEstimator class which utilizes a flux projector to compute a
  smooth flux in an ND space. Provides methods for real valued vectors, and complex vectors.
  The flux is assumed to take the form μ⁻¹∇ × V, where V is either A or E depending on the solver mode.
• Introduces the GradFluxErrorEstimator class which utilizes a flux projector to compute a
  smooth flux in an H1^d space. The flux is assumed to be of the form ϵ ∇ ϕ
• Introduce the CurlFluxCoefficient and GradFluxCoefficient classes which compute the numerical Curl
  and Grad fluxes given a solution.

Issue #, if available:
#2

[sebastiangrimberg]

Looking good @hughcars, I went ahead and made an initial pass as the tests are passing now.

Several style and organizational comments, but one large issue which came up towards the end in the linalg/errorestimator.cpp classes. I see you are allocating a FiniteElementSpace (L2) for the discontinuous flux, and constructing these scalar_mass_matrices and smooth_to_coarse_embed objects to perform the estimation. For reasons noted above, this is something to avoid and I think the loops over elements once the smooth flux GridFunction is computed needs to be reworked. This can be done fully matrix-free, as mfem::GridFunction::ComputeElementLpErrors. There may be good reason you have decided to go your current route vs. the other implementation, but maybe this was done prematurely. I would imagine it looking something like:

... compute mfem::ParGridFunction smooth_flux, and using CurlFluxCoefficient flux_coef ...

for e in 0, ... , NE-1:
    T = mesh.GetElementTransformation(e);
    ir = mfem::IntRules.Get(T.GetGeometryType(), q_order);
    for ip in ir.GetNPoints():
        T.SetIntPoint(ip);
        U = smooth_flux.GetVectorValue(T, ip);
        V = flux_coef.Eval(T, ip);

        // Integrate over the element: \int_e (U-V)^T (U-V) dV = \sum_i w_i ||U_i - V_i||_2^2
        // (w_i = ip.weight * T->Weight(), see GridFunction::ComputeLpError)

        // For the normalization, use \int_e ||U||_2^2 dV = sum_i w_i ||U_i||^2

    end
end

The loop over integration points can use dense matrix arithmetic using the overrides to GetVectorValue and Eval which return a DenseMatrix where each column is the vector at an integration point. I am fairly certain this approach is an improvement in the long run over what you have there, as you avoid the big memory allocations associated with the dense element matrices.

[sebastiangrimberg]

I'm working my way up from the bottom, so have added some (final) comments on ErrorEstimator and ErrorIndicator. I think there is something to be rethought with the ErrorIndicator class which will clarify the normalization and make this class more useful.

I haven't looked at anything inside of drivers/ since I think this changes should be finalized first before I look at the "user" code.

[sebastiangrimberg]

I think we should squash these commits before merging #103 in, and then merging to main. What do you think? If that sounds good feel free to go ahead and do so and rebase and merge #103 here.