Source code for goalie.error_estimation

"""
Tools to automate goal-oriented error estimation.
"""

import firedrake
import ufl
from firedrake import Function, FunctionSpace
from firedrake.functionspaceimpl import WithGeometry
from firedrake.petsc import PETSc

__all__ = ["get_dwr_indicator"]


@PETSc.Log.EventDecorator()
def form2indicator(F):
    r"""
    Given a 0-form, multiply the integrand of each of its integrals by a
    :math:`\mathbb P0` test function and reassemble to give an element-wise error
    indicator.

    Note that a 0-form does not contain any :class:`firedrake.ufl_expr.TestFunction`\s
    or :class:`firedrake.ufl_expr.TrialFunction`\s.

    :arg F: the 0-form
    :type F: :class:`ufl.form.Form`
    :return: the corresponding error indicator field
    :rtype: `firedrake.function.Function`
    """
    if not isinstance(F, ufl.form.Form):
        raise TypeError(f"Expected 'F' to be a Form, not '{type(F)}'.")
    mesh = F.ufl_domain()
    P0 = FunctionSpace(mesh, "DG", 0)
    p0test = firedrake.TestFunction(P0)
    h = ufl.CellVolume(mesh)

    rhs = 0
    for integral in F.integrals_by_type("exterior_facet"):
        ds = firedrake.ds(integral.subdomain_id())
        rhs += h * p0test * integral.integrand() * ds
    for integral in F.integrals_by_type("interior_facet"):
        dS = firedrake.dS(integral.subdomain_id())
        rhs += h("+") * p0test("+") * integral.integrand() * dS
        rhs += h("-") * p0test("-") * integral.integrand() * dS
    for integral in F.integrals_by_type("cell"):
        dx = firedrake.dx(integral.subdomain_id())
        rhs += h * p0test * integral.integrand() * dx

    assert rhs != 0
    indicator = Function(P0)
    firedrake.solve(
        firedrake.TrialFunction(P0) * p0test * firedrake.dx == rhs,
        indicator,
        solver_parameters={
            "snes_type": "ksponly",
            "ksp_type": "preonly",
            "pc_type": "jacobi",
        },
    )
    return indicator


[docs] @PETSc.Log.EventDecorator() def get_dwr_indicator(F, adjoint_error, test_space=None): r""" Given a 1-form and an approximation of the error in the adjoint solution, compute a dual weighted residual (DWR) error indicator. Note that each term of a 1-form contains only one :class:`firedrake.ufl_expr.TestFunction`. The 1-form most commonly corresponds to the variational form of a PDE. If the PDE is linear, it should be written as in the nonlinear case (i.e., with the solution field in place of any :class:`firedrake.ufl_expr.TrialFunction`\s. :arg F: the form :type F: :class:`ufl.form.Form` :arg adjoint_error: a dictionary whose keys are field names and whose values are the approximations to the corresponding components of the adjoint error, or a single such component :type adjoint_error: :class:`firedrake.function.Function` or :class:`dict` with :class:`str` keys and :class:`firedrake.function.Function` values :kwarg test_space: a dictionary whose keys are field names and whose values are the test spaces for the corresponding fields, or a single such test space (or ``None`` to determine the test space(s) automatically) :type test_space: :class:`firedrake.functionspaceimpl.WithGeometry` :returns: the DWR indicator :rtype: :class:`firedrake.function.Function` """ mapping = {} if not isinstance(F, ufl.form.Form): raise TypeError(f"Expected 'F' to be a Form, not '{type(F)}'.") # Process input for adjoint_error as a dictionary if isinstance(adjoint_error, Function): name = adjoint_error.name() if test_space is None: test_space = {name: adjoint_error.function_space()} adjoint_error = {name: adjoint_error} elif not isinstance(adjoint_error, dict): raise TypeError( f"Expected 'adjoint_error' to be a Function or dict, not '{type(adjoint_error)}'." ) # Process input for test_space as a dictionary if test_space is None: test_space = {key: err.function_space() for key, err in adjoint_error.items()} elif isinstance(test_space, WithGeometry): if len(adjoint_error.keys()) != 1: raise ValueError("Inconsistent input for 'adjoint_error' and 'test_space'.") test_space = {key: test_space for key in adjoint_error} elif not isinstance(test_space, dict): raise TypeError( f"Expected 'test_space' to be a FunctionSpace or dict, not '{type(test_space)}'." ) # Construct the mapping for each component for key, err in adjoint_error.items(): if key not in test_space: raise ValueError(f"Key '{key}' does not exist in the test space provided.") fs = test_space[key] if not isinstance(fs, WithGeometry): raise TypeError( f"Expected 'test_space['{key}']' to be a FunctionSpace, not '{type(fs)}'." ) if F.ufl_domain() != err.function_space().mesh(): raise ValueError( "Meshes underlying the form and adjoint error do not match." ) if F.ufl_domain() != fs.mesh(): raise ValueError("Meshes underlying the form and test space do not match.") mapping[firedrake.TestFunction(fs)] = err # Apply the mapping return form2indicator(ufl.replace(F, mapping))