Equations

equations

Generates the UFL form for an equation consisting of individual terms.

This module contains a dataclass to define the structure of mathematical equations within the G-ADOPT library. It provides a convenient way to generate the UFL form required by Firedrake solvers.

Equation(test, trial_space, residual_terms, *, eq_attrs={}, approximation=None, bcs=dict(), quad_degree=None, scaling_factor=1) dataclass

Generates the UFL form for the sum of terms constituting an equation.

The generated UFL form corresponds to a sum of implemented term forms contributing to the equation's residual in the finite element discretisation.

Parameters:

Name	Type	Description	Default
test	Argument | Indexed	Firedrake test function.	required
trial_space	WithGeometry	Firedrake function space of the trial function.	required
residual_terms	InitVar[Callable | list[Callable]]	Equation term or a list thereof contributing to the residual.	required
eq_attrs	InitVar[dict[str, Any]]	Dictionary of fields and parameters used in the equation's weak form.	{}
approximation	BaseApproximation | BaseGIAApproximation | None	G-ADOPT approximation for the system of equations considered.	None
bcs	dict[int, dict[str, Any]]	Dictionary specifying weak boundary conditions (identifier, type, and value).	dict()
quad_degree	InitVar[int | None]	Integer specifying the quadrature degree. If omitted, it is set to 2p + 1, where p is the polynomial degree of the trial space.	None
scaling_factor	Number | Constant	A constant factor used to rescale residual terms.	1

residual(trial)

Generates the UFL form corresponding to the residual terms.

Source code in g-adopt/gadopt/equations.py

def residual(self, trial: fd.Argument | Indexed | fd.Function) -> fd.Form:
    """Generates the UFL form corresponding to the residual terms."""
    return self.scaling_factor * sum(
        term(self, trial) for term in self.residual_terms
    )

cell_edge_integral_ratio(mesh, p)

Ratio C such that \int_f u^2 <= C Area(f)/Volume(e) \int_e u^2 for facets f, elements e, and polynomials u of degree p.

See Equation (3.7), Table 3.1, and Appendix C from Hillewaert's thesis: https://www.researchgate.net/publication/260085826

Source code in g-adopt/gadopt/equations.py

def cell_edge_integral_ratio(mesh: fd.MeshGeometry, p: int) -> int:
    r"""
    Ratio C such that \int_f u^2 <= C Area(f)/Volume(e) \int_e u^2 for facets f,
    elements e, and polynomials u of degree p.

    See Equation (3.7), Table 3.1, and Appendix C from Hillewaert's thesis:
    https://www.researchgate.net/publication/260085826
    """
    match cell_type := mesh.ufl_cell().cellname:
        case "triangle":
            return (p + 1) * (p + 2) / 2.0
        case "quadrilateral" | "interval * interval":
            return (p + 1) ** 2
        case "triangle * interval":
            return (p + 1) ** 2
        case "quadrilateral * interval" | "hexahedron":
            # if e is a wedge and f is a triangle: (p+1)**2
            # if e is a wedge and f is a quad: (p+1)*(p+2)/2
            # here we just return the largest of the the two (for p>=0)
            return (p + 1) ** 2
        case "tetrahedron":
            return (p + 1) * (p + 3) / 3
        case _:
            raise NotImplementedError(f"Unknown cell type in mesh: {cell_type}")

interior_penalty_factor(eq, *, shift=0)

Interior Penalty method

For details on the choice of sigma, see https://www.researchgate.net/publication/260085826

We use Equations (3.20) and (3.23). Instead of getting the maximum over two adjacent cells (+ and -), we just sum (i.e. 2 * avg) and have an extra 0.5 for internal facets.

Source code in g-adopt/gadopt/equations.py

def interior_penalty_factor(eq: Equation, *, shift: int = 0) -> float:
    """Interior Penalty method

    For details on the choice of sigma, see
    https://www.researchgate.net/publication/260085826

    We use Equations (3.20) and (3.23). Instead of getting the maximum over two adjacent
    cells (+ and -), we just sum (i.e. 2 * avg) and have an extra 0.5 for internal
    facets.
    """
    degree = eq.trial_space.ufl_element().degree()
    if not isinstance(degree, int):
        degree = max(degree)

    if degree == 0:  # probably only works for orthogonal quads and hexes
        sigma = 1.0
    else:
        # safety factor: 1.0 is theoretical minimum
        alpha = getattr(eq, "interior_penalty", 2.0)
        num_facets = eq.mesh.ufl_cell().num_facets
        sigma = alpha * cell_edge_integral_ratio(eq.mesh, degree + shift) * num_facets

    return sigma