Antoon's site - Blog

home Blog Flow solver LA and MC Simulations W210 Recent activity

Exact and approximate discrete vector-field operators on a Cartesian mesh

In this blog post, I attempt to introduce exact and accompanying consistent approximate discrete vector-field operators on a Cartesian mesh. These operators will act on and “produce” various styles of discretized fields. So first, the staggering of various field discretizations is introduced.

Field discretizations

We consider a cubic Cartesian mesh of size $L^3$, and origin $(x_0, y_0, z_0)$ that is meshed by cubic cells of size $\Delta^3$. These cells have Cartesian indices $[i,j,k]$. We introduce two discretizations of a scalar field $\phi(x,y,z)$:

• $\phi_v[i,j,k]$ presents the field value at the position of the left-lower-back vertex of the corresponding cell: $$\phi_v[i,j,k] = \phi(x_0+i\Delta, y_0 + j\Delta, z_0 + k\Delta).$$
• $\phi_c[i,j,k]$ presents the field’s cell-averaged value over the corresponding cell: $$\phi_c[i,j,k] = \frac{1}{\Delta^3} \int_{x_0+i\Delta}^{x_0+(i+1)\Delta}\int_{y_0+j\Delta}^{y_0+(j+1)\Delta}\int_{z_0+k\Delta}^{z_0+(k+1)\Delta}\phi(x,y,z)\ \mathrm{d}z\mathrm{d}y\mathrm{d}x.$$

These scalar fields are accompanied by two discretizations for a vector field, $\vec F(x,y,z) = F_x(x,y,z)\vec{e}_x + F_y(x,y,z)\vec{e}_y + F_z(x,y,z)\vec{e}_z$.

• $\vec F_f[i,j,k]$ represents the average of the field’s components over the cell’s face with corresponding normal vector. For example for $F_{x,f}[i,j,k]$, $$F_{x,f}[i,j,k] = \frac{1}{\Delta^2}\int_{y_0+j\Delta}^{y_0+(j+1)\Delta}\int_{z_0+k\Delta}^{z_0+(k+1)\Delta} F_x(x_0 + i\Delta,y,z)\ \mathrm{d}z\mathrm{d}y$$
• $\vec F_e[i,j,k]$ represents the average of the field’s components over the cell’s edge with corresponding tangential vector. For example for $F_{x,e}[i,j,k],$ $$F_{x,e}[i,j,k] = \int_{x_0+i\Delta}^{x_0+(i+1)\Delta} F_x(x,y_0 + j\Delta, z_0 + k\Delta)\mathrm{d}x.$$

Exact $\nabla$ field operators

With these discretizations we can define three exact $\nabla$ operators.

• Exact scalar-gradient operator ($\vec F = \nabla \phi$): We can readily compute the edge-averaged values of $\nabla \phi$ from the vertex values of $\phi$ ($\phi_v$), from the gradient theorem. For example, for $F_x = \frac{\partial \phi}{\partial x}$, $$F_{x,e}[i,j,k] = \frac{1}{\Delta}\left(\phi_v[i + 1, j, k] - \phi_v[i,j,k]\right ).$$ The exact discrete gradient operator is denoted as $\vec F_e = \nabla \phi_v$, where the subscripts, indicating the staggering, also reference the operator definition.
• Exact curl operator ($\vec \omega = \nabla \times \vec F$): We can readily compute the face-averaged curl from an edge-averaged vector field, From Stokes’ theorem. E.g. for $\omega_{x,f}$, $$\omega_{x,k} = \frac{1}{\Delta}\left(F_{y,e}[i,j.k] + F_{z,e}[i,j+1,k] - F_{y,e}[i,j,k+1] - F_{z,e}[i,j,k] \right).$$ The exact discrete curl operator is denoted as $\vec A_f = \nabla \times \vec F_e.$
  • Exact divergence operator ($\phi = \nabla \cdot \vec F$): We can readily compute the cell-averaged divergence ($\nabla \cdot \vec F$) from the face-averaged vector-components, from Gauss’ theorem, $$\phi_c[i,j,k] = \frac{1}{\Delta}\left(F_x[i+1,j,k] + F_y[i,j+1,k] + F_z[i,j,k+1] - F_x[i,j,k] - F_y[i,j,k] + F_z[i,j,k]\right).$$ The exact discrete divergence operator is denoted as $\phi_c = \nabla \cdot \vec F_f$.

Given that these are exact operators, it is no surprise that the two possible “double $\nabla$” identities are satisfied: $$\nabla \times (\nabla \phi_v)_e = \vec 0_f, \ \ \ \text{ and, } \ \ \ \ \ \nabla \cdot (\nabla \times \vec F_e)_f = 0_c.$$ Indeed, if we compute the exact (face-averaged) curl from exact (edge-averaged) gradient (computed from vertex values), we exactly get the zero vector. Further the exact (cell averaged) divergence of a (face-averaged) curl (from an edge-averaged vector field) is also exactly zero. Note that the set of three exact operators always “go up” in dimensionality, $$\text{Vertex-point (0D) field } \xrightarrow{\nabla \phi_v} \text{ Edge-averaged (1D) field} \xrightarrow{\nabla \times \vec F_e} \text{Face-averaged (2D) field } \xrightarrow{\nabla \cdot \vec F_{f}} \text{ Cell-averaged (3D) field}.$$ As such it would be nice if we could augment this set op operator definitions, with operators in the other direction. We call these, now approximate, operators “consistent” when the formulations satisfy the “double $\nabla$” identities exactly.

Consistent approximate $\nabla$ operators.

We introduce three second-order accurate $\nabla$ operators to complement the three exact operators.

• Approximate face-averaged gradient from cell-averaged fields, $$F_{x,f} = \frac{1}{\Delta}\left(\phi_c[i,j,k]-\phi_c[i-1,j,k]\right)$$
• Approximate edge-averaged curl from face-averaged vector fields, $$F_{x,f} = \frac{1}{\Delta}\left(\phi_c[i,j,k]-\phi_c[i-1,j,k]\right)$$
• Approximate vertex-point divergence from edge-averaged fields, $$\phi_{v} = \frac{1}{\Delta}\left( F_{x,e}[i,j,k]-F_{x,e}[i-1,j,k] + F_{y,e}[i,j,k]-F_{y,e}[i-1,j-1,k]+ F_{z,e}[i,j,k]-F_{z,e}[i,j,k - 1]\right)$$

It can be easily verified that these operators satisfy the aforementioned “double $\nabla$” identities exactly, $$\nabla \cdot (\nabla \times \vec F_f)_e = 0 \ \ \ \ \ \ \text{and} \ \ \ \ \ \ \ \nabla \times (\nabla \phi_c)_f = \vec 0,$$ and are therefore labeled “consistent”.

Applications

To finalize this blog post, and to illustrate the use of these formulations, five elementary “applications” for vector-field analysis are introduced.

Application 1: The Laplacian

We can compute the Laplacian ($L = \nabla \cdot \nabla \phi$) of a scalar field “in place”, i.e. from a vertex field, as a vertex field and from a cell-averaged field as a cell-averaged field. Using the exact gradient operator, $$L_v = \nabla \cdot (\nabla \phi_v)_e,$$ or using the exact divergence operator, $$L_c = \nabla \cdot (\nabla \phi_c)_f.$$

Application 2: Two Projection methods

With the Laplacian we can project a vector field onto the space of divergence free vector fields. For a face-averaged vector field ($\vec F_f$) the steps are,

1. Compute the exact cell-averaged divergence: $\nabla \cdot \vec F_f$
2. Solve the discrete Poisson problem for the cell-averaged field $a_c$, $$\nabla \cdot (\nabla a_c)_f = \nabla \cdot \vec F_f$$
3. Reject the gradient of $a_c$ from $F_f$ (containing the divergence) to find the desired divergence-free projected vector field, $\vec P_f$, $$\vec P_f = \vec F_f - \nabla a_c.$$

With this method, one can claim the cell-averaged divergence of $\vec P_f$ is exactly zero. Alternatively, one can project an edge-averaged vectorfield ($\vec F_e$) with similar steps,

1. Compute the approximate vertex-point divergence $\nabla \cdot \vec F_e$
2. Solve the Poisson problem for the vertex-point field $a_v$, $$\nabla \cdot (\nabla a_c)_e = \nabla \cdot \vec F_e$$
3. Reject the gradient of $a_v$ from $F_e$ (containing the divergence) to find the desired projection vector $\vec P_e$, $$\vec P_e = \vec F_e - \nabla a_v.$$

Although the approximate divergence of $\vec P_e$ will be exactly zero, the latter method seems less favorable compared to the former, as for the latter, errors arise for both applications of the approximate divergence operator. I wonder if these two can be related.

Application 3: Vector potential

For a vector field described by a vector potential $\vec F = \nabla \times \vec A$, we can find the vector potential by solving either, $$\nabla \times (\nabla \times \vec A_f)_e = \nabla \times \vec F_e,$$ or, $$\nabla \times (\nabla \times \vec A_e)_f = \nabla \times \vec F_f,$$ depending on the (face or edge) distretization of $\vec F$.

Application 4: Consistent Helmholtz decomposition

Combining the projection (for the scalar potential) and vector potential, we can decompose any scalar field into a (rotation free) scalar gradient and (divergence free) curl field, $$\vec F = \nabla a + \nabla \times \vec A.$$

There are two options:

1. A face-averaged vector field ($\vec F_f$) may be reconstructed by a cell-averaged scalar potential and an edge-averaged vector potential.
2. An edge-averaged vector field ($\vec F_e$) may be reconstructed by a vertex-point scalar potential and a face-averaged vector potential.

Application 5: The vector Laplacian

Since we have not defined a method to compute the gradient of a vector-component field, the vector Laplacian seems illusive. However, using the identity: $$\nabla^2 \vec F = \nabla (\nabla \cdot \vec F) - \nabla\times\nabla\times\vec F,$$ we can approximate the vector Laplacian “in place” for edge-averaged vector fields, $$\vec B_e = \nabla^2 \vec F_e = \nabla (\nabla \cdot \vec F_e)_v - \nabla \times (\nabla \times \vec F_e)_f,$$ and face-averaged vector fields, $$\vec B_f = \nabla^2 \vec F_f = \nabla (\nabla \cdot \vec F_f)_c - \nabla \times (\nabla \times \vec F_f)_e.$$

Further applications of these “applications” will be presented later.

Vector product in $D$ dimensions

It is often¹ said that the vector product (“$\times$”, cross product) cannot be generalized to other than three dimensions. But it depends a bit on how one defines the vector product. If we forgo the requirement that is must be a binary operation (i.e. the product of two vectors) but instead require $D - 1$ vectors for the $D$ dimensional vector product, it generalizes quite naturally: $$\vec u = \times(\vec v_1, \vec v_2, ... , \vec v_{D-1}) = \vec{v}_1\times\vec v_2 \times ..., \vec v_{D-1} = \begin{vmatrix} \vec e_1 & \vec e_2 & ... & \vec e_D\\ v_{1,1} & \vec v_{1,2} & ... & v_{1,D}\\ v_{2,1} & \vec v_{2,2} & ... & v_{2,D}\\ \vdots & \vdots & \ddots & \vdots\\ v_{D-1,1} & \vec v_{D-1,2} & ... & v_{D-1,D}\end{vmatrix}.$$ By the properties of determinants, it remains to be distributive over addition (i.e. it is a linear, product-style operator), $\vec u = \vec 0$ if any two vectors in the product are parallel, and $\vec u$ is perpendicular to every argument of the product. Using this definition, we can use it for any dimension larger than one. For example, the vector product of $\vec V_1 = 2\vec e_1 + 1\vec e_2$ in 2 dimensions is $\vec U$, $$\vec U = \times(\vec V_1) = \begin{vmatrix}\vec e_1 & \vec e_2\\ 2& 1\end{vmatrix} = 1\vec{e}_1 - 2\vec{e}_2.$$

---

1. It is not often a topic of discussion, really.↩︎