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 Ex1-Volume, 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 Area 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

(8)\[ \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

(9)\[ \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 (8), 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 Area 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. 6). This problem can be run with:

./area -problem sphere

Fig. 6 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. 7), to \(\bm{X}=(X,Y) \in \textrm{I}\) on the reference element, via the chain rule

(10)\[ \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

(11)\[ \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. 7 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 (10), 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 bakeoff problems (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 (11), are the same as in the Sphere example. For the Poisson’s problem, BP3-BP6, on the cubed-sphere, in addition to equation (11), the pseudo-inverse of \(\partial \overset{\circ}{\bm{x}} / \partial \bm{X}\) is used to derive the contravariant metric tensor (please see figure Fig. 7 for a reference of the notation used). We begin by expressing the Moore-Penrose (left) pseudo-inverse:

(12)\[ \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

(13)\[ \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 (13) 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 (14), using an unstructured high-order finite element discretization. All of the operators associated with the geometric multigrid are implemented in libCEED.

(14)\[ -\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 Operator Decomposition, 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.