# PETSc demos and BPs¶

## Area¶

This example is located in the subdirectory examples/petsc. It demonstrates a simple usage of libCEED with PETSc to calculate the surface area of a closed surface. The code uses higher level communication protocols for mesh handling in PETSc’s DMPlex. This example has the same mathematical formulation as , with the exception that the physical coordinates for this problem are $$\bm{x}=(x,y,z)\in \mathbb{R}^3$$, while the coordinates of the reference element are $$\bm{X}=(X,Y) \equiv (X_0,X_1) \in \textrm{I} =[-1,1]^2$$.

### Cube¶

This is one of the test cases of the computation of the of a 2D manifold embedded in 3D. This problem can be run with:

./area -problem cube


This example uses the following coordinate transformations for the computation of the geometric factors: from the physical coordinates on the cube, denoted by $$\bar{\bm{x}}=(\bar{x},\bar{y},\bar{z})$$, and physical coordinates on the discrete surface, denoted by $$\bm{{x}}=(x,y)$$, to $$\bm{X}=(X,Y) \in \textrm{I}$$ on the reference element, via the chain rule

(5)$\frac{\partial \bm{x}}{\partial \bm{X}}_{(2\times2)} = \frac{\partial {\bm{x}}}{\partial \bar{\bm{x}}}_{(2\times3)} \frac{\partial \bar{\bm{x}}}{\partial \bm{X}}_{(3\times2)},$

with Jacobian determinant given by

(6)$\left| J \right| = \left\|col_1\left(\frac{\partial \bar{\bm{x}}}{\partial \bm{X}}\right)\right\| \left\|col_2 \left(\frac{\partial \bar{\bm{x}}}{\partial \bm{X}}\right) \right\|$

We note that in equation (5), the right-most Jacobian matrix $${\partial\bar{\bm{x}}}/{\partial \bm{X}}_{(3\times2)}$$ is provided by the library, while $${\partial{\bm{x}}}/{\partial \bar{ \bm{x}}}_{(2\times3)}$$ is provided by the user as

$\left[ col_1\left(\frac{\partial\bar{\bm{x}}}{\partial \bm{X}}\right) / \left\| col_1\left(\frac{\partial\bar{\bm{x}}}{\partial \bm{X}}\right)\right\| , col_2\left(\frac{\partial\bar{\bm{x}}}{\partial \bm{X}}\right) / \left\| col_2\left(\frac{\partial\bar{\bm{x}}}{\partial \bm{X}}\right)\right\| \right]^T_{(2\times 3)}.$

### Sphere¶

This problem computes the surface of a tensor-product discrete sphere, obtained by projecting a cube inscribed in a sphere onto the surface of the sphere. This discrete surface is sometimes referred to as a cubed-sphere (an example of such as a surface is given in figure Fig. 4). This problem can be run with:

./area -problem sphere


Fig. 4 Example of a cubed-sphere, i.e., a tensor-product discrete sphere, obtained by projecting a cube inscribed in a sphere onto the surface of the sphere.

This example uses the following coordinate transformations for the computation of the geometric factors: from the physical coordinates on the sphere, denoted by $$\overset{\circ}{\bm{x}}=(\overset{\circ}{x},\overset{\circ}{y},\overset{\circ}{z})$$, and physical coordinates on the discrete surface, denoted by $$\bm{{x}}=(x,y,z)$$ (depicted, for simplicity, as coordinates on a circle and 1D linear element in figure Fig. 5), to $$\bm{X}=(X,Y) \in \textrm{I}$$ on the reference element, via the chain rule

(7)$\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}_{(3\times2)} = \frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{x}}_{(3\times3)} \frac{\partial\bm{x}}{\partial \bm{X}}_{(3\times2)} ,$

with Jacobian determinant given by

(8)$\left| J \right| = \left| col_1\left(\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}\right) \times col_2 \left(\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}\right)\right| .$

Fig. 5 Sketch of coordinates mapping between a 1D linear element and a circle. In the case of a linear element the two nodes, $$p_0$$ and $$p_1$$, marked by red crosses, coincide with the endpoints of the element. Two quadrature points, $$q_0$$ and $$q_1$$, marked by blue dots, with physical coordinates denoted by $$\bm x(\bm X)$$, are mapped to their corresponding radial projections on the circle, which have coordinates $$\overset{\circ}{\bm{x}}(\bm x)$$.

We note that in equation (7), the right-most Jacobian matrix $${\partial\bm{x}}/{\partial \bm{X}}_{(3\times2)}$$ is provided by the library, while $${\partial \overset{\circ}{\bm{x}}}/{\partial \bm{x}}_{(3\times3)}$$ is provided by the user with analytical derivatives. In particular, for a sphere of radius 1, we have

$\overset{\circ}{\bm x}(\bm x) = \frac{1}{\lVert \bm x \rVert} \bm x_{(3\times 1)}$

and thus

$\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{x}} = \frac{1}{\lVert \bm x \rVert} \bm I_{(3\times 3)} - \frac{1}{\lVert \bm x \rVert^3} (\bm x \bm x^T)_{(3\times 3)} .$

## Bakeoff problems and generalizations¶

The PETSc examples in this directory include a full suite of parallel (BPs) using a “raw” parallel decomposition (see bpsraw.c) and using PETSc’s DMPlex for unstructured grid management (see bps.c). A generalization of these BPs to the surface of the cubed-sphere are available in bpssphere.c.

### Bakeoff problems on the cubed-sphere¶

For the $$L^2$$ projection problems, BP1-BP2, that use the mass operator, the coordinate transformations and the corresponding Jacobian determinant, equation (8), are the same as in the example. For the Poisson’s problem, BP3-BP6, on the cubed-sphere, in addition to equation (8), the pseudo-inverse of $$\partial \overset{\circ}{\bm{x}} / \partial \bm{X}$$ is used to derive the contravariant metric tensor (please see figure Fig. 5 for a reference of the notation used). We begin by expressing the Moore-Penrose (left) pseudo-inverse:

(9)$\frac{\partial \bm{X}}{\partial \overset{\circ}{\bm{x}}}_{(2\times 3)} \equiv \left(\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}\right)_{(2\times 3)}^{+} = \left(\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}_{(2\times3)}^T \frac{\partial\overset{\circ}{\bm{x}}}{\partial \bm{X}}_{(3\times2)} \right)^{-1} \frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}_{(2\times3)}^T .$

This enables computation of gradients of an arbitrary function $$u(\overset{\circ}{\bm x})$$ in the embedding space as

$\frac{\partial u}{\partial \overset{\circ}{\bm x}}_{(1\times 3)} = \frac{\partial u}{\partial \bm X}_{(1\times 2)} \frac{\partial \bm X}{\partial \overset{\circ}{\bm x}}_{(2\times 3)}$

and thus the weak Laplacian may be expressed as

(10)$\int_{\Omega} \frac{\partial v}{\partial \overset\circ{\bm x}} \left( \frac{\partial u}{\partial \overset\circ{\bm x}} \right)^T \, dS = \int_{\Omega} \frac{\partial v}{\partial \bm X} \underbrace{\frac{\partial \bm X}{\partial \overset\circ{\bm x}} \left( \frac{\partial \bm X}{\partial \overset\circ{\bm x}} \right)^T}_{\bm g_{(2\times 2)}} \left(\frac{\partial u}{\partial \bm X} \right)^T \, dS$

where we have identified the $$2\times 2$$ contravariant metric tensor $$\bm g$$ (sometimes written $$\bm g^{ij}$$), and where now $$\Omega$$ represents the surface of the sphere, which is a two-dimensional closed surface embedded in the three-dimensional Euclidean space $$\mathbb{R}^3$$. This expression can be simplified to avoid the explicit Moore-Penrose pseudo-inverse,

$\bm g = \left(\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}^T \frac{\partial\overset{\circ}{\bm{x}}}{\partial \bm{X}} \right)^{-1}_{(2\times 2)} \frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}_{(2\times3)}^T \frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}_{(3\times2)} \left(\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}^T \frac{\partial\overset{\circ}{\bm{x}}}{\partial \bm{X}} \right)^{-T}_{(2\times 2)} = \left(\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}^T \frac{\partial\overset{\circ}{\bm{x}}}{\partial \bm{X}} \right)^{-1}_{(2\times 2)}$

where we have dropped the transpose due to symmetry. This allows us to simplify (10) as

$\int_{\Omega} \frac{\partial v}{\partial \overset\circ{\bm x}} \left( \frac{\partial u}{\partial \overset\circ{\bm x}} \right)^T \, dS = \int_{\Omega} \frac{\partial v}{\partial \bm X} \underbrace{\left(\frac{\partial \overset{\circ}{\bm{x}}}{\partial \bm{X}}^T \frac{\partial\overset{\circ}{\bm{x}}}{\partial \bm{X}} \right)^{-1}}_{\bm g_{(2\times 2)}} \left(\frac{\partial u}{\partial \bm X} \right)^T \, dS ,$

which is the form implemented in qfunctions/bps/bp3sphere.h.

## Multigrid¶

This example is located in the subdirectory examples/petsc. It investigates $$p$$-multigrid for the Poisson problem, equation (11), using an unstructured high-order finite element discretization. All of the operators associated with the geometric multigrid are implemented in libCEED.

(11)$-\nabla\cdot \left( \kappa \left( x \right) \nabla x \right) = g \left( x \right)$

The Poisson operator can be specified with the decomposition given by the equation in figure , and the restriction and prolongation operators given by interpolation basis operations, $$\bm{B}$$, and $$\bm{B}^T$$, respectively, act on the different grid levels with corresponding element restrictions, $$\bm{G}$$. These three operations can be exploited by existing matrix-free multigrid software and smoothers. Preconditioning based on the libCEED finite element operator decomposition is an ongoing area of research.