Helmholtz Decomposition of Vector Fields Using an Optimal Preconditioned Conjugate Gradient Algorithm

Jorge Lopez

UJAT, Tabasco, México.

DOI: 10.4236/jamp.2023.115086                

Abstract

In this article, we study numerically a Helmholtz decomposition methodology, based on a formulation of the mathematical model as a saddle-point problem. We use a preconditioned conjugate gradient algorithm, applied to an associated operator equation of elliptic type, to solve the problem. To solve the elliptic partial differential equations, we use a second order mixed finite element approximation for discretization. We show, using 2-D synthetic vector fields, that this approach, yields very accurate solutions at a low computational cost compared to traditional methods with the same order of approximation.

Keywords

Helmholtz Decomposition, Self-Adjoint Operator, Optimal Preconditioning, Finite Element

1. Introduction

Several problems and applications require a deep knowledge of a vector field over a region. An important tool in this context is the Helmholtz decomposition theorem (HDT) that guarantees, under certain hipothesis [1] [2] , a unique decomposition of a vector field $u$ defined in a domain $\Omega$ into its solenoidal $u_s$ and irrotational $u_p$ components. This is,

$u=u_s+u_p\mathrm{,}\text{in}\Omega \mathrm{.}$ (1)

This decomposition is important, for example, because the properties like incompressibility and vorticity can be studied directly by studying these components.

The application areas of the HD include (but are not restricted to) electromagnetism, linear elasticity, fluid mechanic [3] , astrophysics [4] , geophysics [5] , computer vision and robotic [6] [7] [8] , image processesing [9] . The HDT also establishes that the Helmoltz decomposition components can be obtained in terms of a vector potential (the solenoidal) and a scalar potential (the irrotational), solution of certain Poisson equations with different kinds of boundary conditions, one of them being $u_s\cdot n=0$ , on $\Gamma =\partial \Omega$ . These elliptic equations cannot be solved explicitly, so it is necessary to apply numerical methods to solve it. The most common mesh dependant methods are the finite element methods (FEM) [6] [10] [11] , Finite difference methods (FDM) [9] [12] , Galerkin Formulation Using Finite Differences [13] , Fourier and Wavelets Domains [14] [15] [16] . Among the meshsless methods, we can mention the Smoothed Particle Hydrodynamics (SPH) [17] [18] [19] and the Statistical Learning Using Matrix-Valued Kernels [20] [21] [22] .

The problem we are considering in this work is the numerical decomposition, in a given domain $\Omega$ , of a vector field $u$ into its solenoidal and irrotational components with the more general boundary condition: $u_s\cdot n=0$ is satisfied only in ${\Gamma }_N$ , a section of $\Gamma$ where the solenoidal vector field is tangential to the boundary. This kind of boundary conditions appears, for example when we want to adjust vector fields from experimental or incomplete data. Let u be given and defined in the Lipschitz bounded domain $\Omega$ , with boundary $\Gamma ={\Gamma }_N\cup {\Gamma }_D$ . We want to solve numerically the following model:

$u=u_s+u_p\mathrm{,}\text{in}\Omega \mathrm{,}$ (2)

$\nabla \cdot u_s=\mathrm{0,}\text{in}\Omega \mathrm{,}$ (3)

$u_p=\nabla p\text{in}\Omega \mathrm{,}$ (4)

$u_s\cdot n=\mathrm{0,}\text{on}{\Gamma }_N\mathrm{,}$ (5)

where p is a scalar potential, $n$ is the unitary outer vector, normal to the boundary $\Gamma$ and ${\Gamma }_D=\Gamma -{\Gamma }_N$ .

We apply a methodology based on the solution of saddle point problems to find simultaneously the pair $\left(u_s\mathrm{,}p\right)$ (and consequently $u_p=u-u_s$ ) by using an optimal preconditioned conjugate gradient algorithm combined with a mixed finite element method to solve elliptic problems. The associated stiffness matrix components are aproximated by a trapesoidal rule that produce a diagonal matrix so that the solution of the stiffness system is easy. This property together with the use of an optimal preconditioned resulted in a very fast algorithm.

The structure of this article is as follows: In Section 2, we consider the mathematical formulation of the problem and we describe the classical elliptic problems for scalar potentials allowing the decomposition. In Section 3, we describe the variational version of the problem. After reformulating the problem, we describe a preconditioned conjugate gradient (PCG-algorithm), where a mixed finite element method is used to solve the elliptic subproblems at each iteration. In Section 4, we present and discuss the numerical results, and finally, in Section 5, we give some concluding remarks.

2. Mathematical Formulation of the Helmholtz Decomposition

The Helmholtz decomposition theorem [2] , establishes that if $\Omega \subset {ℝ}^d\left(d\ge 2\right)$ is a bounded Lipschitz domain with boundary $\Gamma$ , then, each $u\in L_2\left(\Omega \right)$ have a unique Helmholtz decomposition

$u=u_s+u_p\mathrm{,}$ (6)

if we know that the solenoidal component satisfies $u_s\cdot n=0$ in $\Gamma$ , i.e. $u_s\in H_{div}\left(\Omega \right)$ and $u_p\in G\left(\Omega \right)$ , where

$H_{div}\left(\Omega \right)=\left\{w\in H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)\mathrm{<mo>\; :\; </mo>}\nabla \cdot w=\mathrm{0\; <mo>\; \&\; </mo>}w\cdot n=0\text{in}\Gamma \right\}\mathrm{,}$ (7)

$H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)=\left\{v\in L_2\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}\nabla \cdot v\in L_2\left(\Omega \right)\right\}$ (8)

and

$G\left(\Omega \right)=\left\{v\in L_2\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}\exists p\in L_2\left(\Omega \right)\text{such}\text{that}v=\nabla p\right\}\mathrm{.}$ (9)

Spaces $H_{div}\left(\Omega \right)$ and $G\left(\Omega \right)$ are orthogonal in $L_2\left(\Omega \right)$ .

If $u\in H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)$ and the boundary condition $u_s\cdot n=0$ in $\Gamma$ is given, the irrotational vector field $u_p=\nabla p$ can be calculated by obtaining the scalar field $p\left(x\right)\in H^1\left(\Omega \right)\mathrm{<mo>\; /\; </mo>}ℝ$ , as solution of the Poisson problem with Neumann boundary condition

$\begin{array}{l}-\Delta p=-\nabla \cdot u\mathrm{,}\text{in}\Omega \mathrm{,}\\ \nabla p\cdot n=u\cdot n\mathrm{,}\text{on}\Gamma \mathrm{,}\end{array}$ (10)

where $H^1\left(\Omega \right)\mathrm{<mo>\; /\; </mo>}ℝ=\left\{v\in H^1\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}{\displaystyle {\int }_{\Omega }}v\left(x\right)\text{d}x=0\right\}$ .

If we first solve the elliptic problem (10) to find $p\left(x\right)$ then we calculate $u_p=\nabla p$ and finally we calculate $u_s=u-u_p$ .

In order to describe our methodology for finding the Helmholtz decomposition, we consider a given vector field $u\in H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)$ satisfying (2)-(5). From (2), we have

$u-u_s=u_p\mathrm{,}$

or

$u-u_s=\nabla p\mathrm{,}\text{for}\text{some}p\in L_2\left(\Omega \right)\mathrm{,}$

From this,

${\displaystyle {\int }_{\Omega }}\left(u-u_s\right)\cdot v\text{d}x={\displaystyle {\int }_{\Omega }}\nabla p\cdot v\text{d}x\mathrm{,}\forall v\in H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)\mathrm{.}$

Since

${\displaystyle {\int }_{\Omega }}\nabla \lambda \cdot v\text{d}x={\displaystyle {\int }_{\Gamma }}\lambda v\cdot n\text{d}\gamma -{\displaystyle {\int }_{\Omega }}\lambda \nabla \cdot v\text{d}x\mathrm{,}\forall v\in H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)\mathrm{,}\forall \lambda \in L_2\left(\Omega \right)\mathrm{,}$

we can write (if $p=0$ in ${\Gamma }_D$ and $v\cdot n=0$ on ${\Gamma }_N$ )

${\displaystyle {\int }_{\Omega }}\left(u-u_s\right)\cdot v\text{d}x=-{\displaystyle {\int }_{\Omega }}p\nabla \cdot v\text{d}x\mathrm{,}\forall v\in V_N\mathrm{,}$

where $V_N$ is the space of vector fields defined as

$V_N=\left\{v\in H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)\mathrm{<mo>\; :\; </mo>}v\cdot n=0\text{on}{\Gamma }_N\right\}\mathrm{.}$ (11)

Considering condition (3) we conclude that the pair $\left(u_s\mathrm{,}p\right)\in V_N\times L_2\left(\Omega \right)$ , satisfies

$\{\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left(u-u_s\right)\cdot v\text{d}x=-{\displaystyle {\int }_{\Omega }}p\nabla \cdot v\text{d}x\mathrm{,}\forall v\in V_N\mathrm{,}\\ {\displaystyle {\int }_{\Omega }}q\nabla \cdot u_s\text{d}x=\mathrm{0,}\forall q\in L_2\left(\Omega \right)\mathrm{.}\end{array}$ (12)

Next section is devoted to find the pair $\left(u_s\mathrm{,}p\right)\in V_N\times L_2\left(\Omega \right)$ , satisfying (12).

3. A Preconditioned Conjugate Gradient Algorithm

Following [23] , assume that $\left(u_s\mathrm{,}p\right)$ is solution of problem (12) with

$u-u_s=u_p\mathrm{.}$ (13)

Considering that $\nabla \cdot u_s=0$ , it follows that $\nabla \cdot u_p=\nabla \cdot u$ , so we have that p satisfies the following operational equation

$Ap=\nabla \cdot u\mathrm{,}$ (14)

where $A\mathrm{<mo>\; :\; </mo>}{L}_2\left(\Omega \right)\to L_2\left(\Omega \right)$ is the operator defined by

$A\mu =\nabla \cdot u_{\mu }\mathrm{,}$ (15)

with $u_{\mu }$ the solution of

$\{\begin{array}{l}{u}_{\mu }\in H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)\mathrm{,}\\ {\displaystyle {\int }_{\Omega }}{u}_{\mu }\cdot v\text{d}x=-{\displaystyle {\int }_{\Omega }}\mu \nabla \cdot v\text{d}x\mathrm{,}\forall v\in V_N\mathrm{.}\end{array}$ (16)

Note that subscript s in (13) have not the same meaning as p and $\mu$ in (16). The properties of operator A (linear, self-adjoint and strongly elliptic, [23] ), allow us to solve equation (14) by using a conjugate gradient algorithm operating in $L_2\left(\Omega \right)$ . Also, we can use the optimal preconditioner $B^{-1}\mathrm{<mo>\; :\; </mo>}{L}_2\left(\Omega \right)\to L_2\left(\Omega \right)$ defined by

$Bq={\varphi }_q,$ (17)

where ${\varphi }_q$ solves the problem

${\displaystyle {\int }_{\Omega }}\nabla {\varphi }_q\cdot \nabla \psi \text{d}x={\displaystyle {\int }_{\Omega }}q\psi \text{d}x\mathrm{,}\forall \psi \in H_D^1\left(\Omega \right)\mathrm{,}$ (18)

${\varphi }_q=0\text{on}{\Gamma }_D\mathrm{,}$ (19)

$\nabla {\varphi }_q\cdot n=0\text{on}{\Gamma }_N\mathrm{,}$ (20)

with

$H_D^1\left(\Omega \right)=\left\{\psi \in H^1\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}\psi =0\text{on}{\Gamma }_D\right\}\mathrm{.}$ (21)

Next we enumerate the three stages of the preconditioned conjugate gradient algorithm to solve the equation $Ap=b$ :

1) Initialization: $p^0$ given, $g^0=A{p}^0-b$ , ${\hat{g}}^0=B{g}^0$ , $d^0=-{\hat{g}}^0$ .

2) Descent: For $m\ge 0$ , assuming we know $p^m\mathrm{,}{g}^m\mathrm{,}{\hat{g}}^m\mathrm{,}{d}^m$ , find $p^{m+1}\mathrm{,}{g}^{m+1},{\hat{g}}^{m+1}\mathrm{,}{d}^{m+1}$ by

$p^{m+1}=p^m+{\alpha }_m{d}^m$ where ${\alpha }_m=\langle g^m\mathrm{,}{\hat{g}}^m\rangle /\langle d^m\mathrm{,}A{d}^m\rangle$ .

$g^{m+1}=g^m+{\alpha }_{m}A{d}^m\mathrm{,}$

${\hat{g}}^{m+1}={\hat{g}}^m+{\alpha }_{m}B\left(A{d}^m\right)\mathrm{.}$

3) Test of convergence and new conjugate direction:

If $\langle g^{m+1}\mathrm{,}{\hat{g}}^{m+1}\rangle \le \epsilon \langle g^0\mathrm{,}{\hat{g}}^0\rangle$ , take $p=p^{m+1}$ and stop.

Else $d^{m+1}=-{\hat{g}}^{m+1}+{\beta }_m{d}^m$ with ${\beta }_m=\frac{\langle g^{m+1}\mathrm{,}{\hat{g}}^{m+1}\rangle }{\langle g^m\mathrm{,}{\hat{g}}^m\rangle }$

Do $m=m+1$ and return to 2.

Using Equations (15)-(16), which define operator A, and Equations (17)-(18) which define operator B, the detailed preconditioned conjugate gradient algorithm (PCG) to obtain p and $u_s$ (and $u_p$ ) is as follows:

Initialization

1) Given $p^0\in L_2\left(\Omega \right)$ , solve

$\{\begin{array}{l}{u}_p^0\in H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)\mathrm{,}\\ {\displaystyle {\int }_{\Omega }}\left(u_p^0\right)\cdot v\text{d}x=-{\displaystyle {\int }_{\Omega }}{p}^0\nabla \cdot v\text{d}x\mathrm{,}\forall v\in V_N\mathrm{.}\end{array}$

2) Let $g^0=-\nabla \cdot u^0$ , where $u^0=-\left(u_p^0-u\right)$ .

3) Solve

$\{\begin{array}{l}{\varphi }^0\in H_D^1\left(\Omega \right)\mathrm{,}\\ {\displaystyle {\int }_{\Omega }}\left(\nabla {\varphi }^0\right)\cdot \nabla \psi \text{d}x={\displaystyle {\int }_{\Omega }}{g}^0\psi \text{d}x\mathrm{,}\forall \psi \in H_D^1\left(\Omega \right)\mathrm{.}\end{array}$

4) Let ${\hat{g}}^0={\varphi }^0$ , $d^0=-{\hat{g}}^0$ .

Descent

For $m\ge 0$ , assuming $p^m$ , $g^m$ , ${\hat{g}}^m$ , $d^m$ , $u^m$ are known, compute $p^{m+1}$ , $g^{m+1}$ , ${\hat{g}}^{m+1}$ , $d^{m+1}$ and $u^{m+1}$ , using the following steps:

5) Solve

$\{\begin{array}{l}{u}^m\in H\left(div\mathrm{<mo>\; ;\; </mo>}\Omega \right)\mathrm{,}\\ {\displaystyle {\int }_{\Omega }}\left(u^m\right)\cdot v\text{d}x=-{\displaystyle {\int }_{\Omega }}{\text{d}}^m\nabla \cdot v\text{d}x\mathrm{,}\forall v\in V_N\mathrm{.}\end{array}$

6) Let ${\bar{g}}^m=\nabla \cdot u^m$ .

7) Let ${\alpha }_m={\displaystyle {\int }_{\Omega }}{g}^m{\hat{g}}^m\text{d}x/{\displaystyle {\int }_{\Omega }}{\bar{g}}^m{d}^m\text{d}x$ .

8) Solve

$\{\begin{array}{l}{\bar{\varphi }}^m\in H_D^1\left(\Omega \right)\mathrm{,}\\ {\displaystyle {\int }_{\Omega }}\left(\nabla {\bar{\varphi }}^m\right)\cdot \nabla \psi \text{d}x={\displaystyle {\int }_{\Omega }}{\bar{g}}^m\psi \text{d}x\mathrm{,}\forall \psi \in H_D^1\left(\Omega \right)\mathrm{.}\end{array}$

9) Set

$p^{m+1}=p^m+{\alpha }_m{d}^m\mathrm{,}$

$u^{m+1}=u^m-{\alpha }_m{u}^m\mathrm{,}$

$g^{m+1}=g^m+{\alpha }_m{\bar{g}}^m\mathrm{,}$

${\hat{g}}^{m+1}={\hat{g}}^m+{\alpha }_m{\bar{\varphi }}^m\mathrm{.}$

Test of convergence and new descent direction

If ${\displaystyle {\int }_{\Omega }}{g}^{m+1}{\hat{g}}^{m+1}\text{d}x/{\displaystyle {\int }_{\Omega }}{g}^0{\hat{g}}^0\text{d}x<\epsilon$ , then do $p=p^{m+1}$ , $u_s=u^{m+1}$ and stop. Otherwise, do the following:

10) Compute ${\beta }_m={\displaystyle {\int }_{\Omega }}{g}^{m+1}{\hat{g}}^{m+1}\text{d}x/{\displaystyle {\int }_{\Omega }}{g}^m{\hat{g}}^m\text{d}x$

11) Set $d^{m+1}=-{\hat{g}}^{m+1}+{\beta }_m{d}^m$ .

12) Do $m=m+1$ and return to 5.

Note that in this algorithm, $u_s$ and p are computed simultaneously.

To approximate the functions belonging to the spaces $V_N$ and $L_2\left(\Omega \right)$ , we make use of the Bercovier-Pironneau finite element approximation using two triangulation of $\Omega$ , the coarse mesh ${\mathcal{T}}_{2h}$ and the fine mesh ${\mathcal{T}}_h$ obtained from the coarse triangulation through a regular subdivision of each triangle $T\in {\mathcal{T}}_{2h}$ , as shown in Figure 1.

The vector-valued functions such as $u_s$ are approximated in the fine mesh and the scalar functions, such as p are approximated in the coarse mesh. Since $u_s$ is obtained on the fine mesh, its resolution is the same as that obtained when solving the elliptics problems (10) in ${\mathcal{T}}_h$ .

To measure the global difference between the exact solenoidal field $u_s$ and the computed solenoidal field $u_{sh}$ we take the relative error

$e_{sr}=\frac{{\Vert u_s-u_{sh}\Vert }_2}{{\Vert u_s\Vert }_2}\mathrm{.}$ (22)

Also, we computed the $L^2$ -norm of the divergence of $u_{sh}$ , ${\Vert \nabla \cdot u_{sh}\Vert }_{L_2\left(\Omega \right)}$ , which we denote by ndiv in the next section.

4. Numerical Results

To show the performance of the PCG-algorithm we chose two synthetic vector fields. In Example 1 we consider a vector field satisfying the boundary condition

$u_s\cdot n=0\text{,}\text{on}\Gamma \mathrm{.}$ (23)

In Example 2 we consider a vector field satisfying the boundary condition

$u_s\cdot n=\mathrm{0,}$ (24)

but only in ${\Gamma }_N$ , a section of $\Gamma$ .

Figure 1. Element in ${\mathcal{T}}_{2h}$ : triangle ABC. Elements in ${\mathcal{T}}_h$ : triangles AQP, PRC, PQR and QBR.

All numerical calculations were done in a Toshiba PC: Portege R705, Windows 7, 64 Bits, Intel Processor CORE i3, 2.27 GHz, 3 GB Ram.

Example 1. We consider the 2D vector field

$u\left(x\mathrm{,}y\right)=u_s\left(x\mathrm{,}y\right)+u_c\left(x\mathrm{,}y\right)\mathrm{,}\Omega =\left(-\mathrm{1,1}\right)\times \left(-\mathrm{1,1}\right)\mathrm{,}$ (25)

for which the solenoidal and irrotational components are

$u_s\left(x\mathrm{,}y\right)=\left(-\sin \left(\pi x\right)\cos \left(\pi y\right),\cos \left(\pi x\right)\sin \left(\pi y\right)\right),$ (26)

$u_c\left(x\mathrm{,}y\right)=\left(\pi \cos \left(\pi \left(x+y\right)\right),\pi \cos \left(\pi \left(x+y\right)\right)\right)=\nabla \left(\sin \left(\pi \left(x+y\right)\right)\right)\mathrm{.}$ (27)

In this case,

$u_s\cdot n=0\text{in}{\Gamma }_N=\Gamma \mathrm{.}$ (28)

Consequently, ${\Gamma }_D=\varnothing$ . To solve the elliptic problems in the PCG algorithm we use simultaneously a coarse triangular mesh ${\mathcal{T}}_{2h}$ (blue in Figure 2) to approximate the scalar potential p and a fine triangular mesh ${\mathcal{T}}_{2h}$ (red in Figure 2) to approximate the solenoidal component.

In Figure 3, we show the exact vector fields $u\left(x\mathrm{,}y\right)$ , $u_s\left(x\mathrm{,}y\right)$ and $u_c\left(x\mathrm{,}y\right)$ .

We want to see how well we can approximate the solenoidal component and the scalar potential p applying the PCG algorithm. A summary of the numerical results is shown in Table 1, where we show the relative error $e_{sr}$ of the solenoidal component, the $L^2$ -norm of the divergence ndiv of the solenoidal component, as well as the number of iterations to get convergence, up to the given tolerance ( $\epsilon ={10}^{-4}$ ), for several mesh sizes.

Figure 2. An example of a 5 × 5 coarse mesh and a 9 × 9 refined mesh. The scalar potential p is approximated in the coarse mesh ${\mathcal{T}}_{2h}$ (blue). The solenoidal component $u_s$ is approximated in the fine mesh ${\mathcal{T}}_h$ (red). On each blue triangle there are four red triangles.

Figure 3. Exact vector fields for Example 1. Left: $u$ . Middle: solenoidal component $u_s$ . Right: irrotational component $u_c$ .

Table 1. Numerical results for Example 1. $\epsilon ={10}^{-4}$ .

In Figure 4, we show the stream lines for the exact (left) and aproximated (right) solenoidal vector fields. In Figure 5, we show the exact (left) and approximated (right) p.

It is important to mention that the error bounds in Table 1 were obtained using in the stop criterion of the conjugate gradient algorithm the relaxed value $\epsilon ={10}^{-4}$ and the results can be improved if we use a more restrictive $\epsilon$ , as is shown in the next example.

Example 2. Here we consider the 2D vector field $u\left(x\mathrm{,}y\right)=(\; x\; ,0\; )$

$u\left(x\mathrm{,}y\right)=u_s\left(x\mathrm{,}y\right)+u_c\left(x\mathrm{,}y\right)=\left(x\mathrm{,0}\right)\mathrm{,}\Omega =\left(\mathrm{1,2}\right)\times \left(\mathrm{0,1}\right)\mathrm{,}$ (29)

for which the solenoidal and irrotational components are

$u_s\left(x\mathrm{,}y\right)=\left(x\mathrm{,}-y\right)\mathrm{,}$ (30)

$u_c\left(x\mathrm{,}y\right)=\left(\mathrm{0,}y\right)=\nabla p\left(x\mathrm{,}y\right)=\nabla \left(y^2/2+c\right)\mathrm{,}c\text{any}\text{constant}\mathrm{.}$ (31)

In this example,

$u_s\cdot n=0\text{on}{\Gamma }_N=\left\{\left(x\mathrm{,}y\right)\in \partial \Omega \mathrm{<mo>\; :\; </mo>}y=0\right\}\mathrm{,}$ (32)

so that, ${\Gamma }_D$ is the top and vertical boundary of $\Omega$ . Our formulation implies that $p=0$ in ${\Gamma }_D$ , but this is only possible in the top part of ${\Gamma }_D$ , and we do not want (and we are not able) to impose additional boundary conditions in p, so we impose exact boundary conditions in $u_s$ in the vertical boundary. Table 2 shows the numerical results obtained for this test problem.

For comparison reasons we have repeated the calculations in example 2 but using a more restrictive stop tolerance of $\epsilon ={10}^{-12}$ for the conjugate gradient

Figure 4. Stream lines for the exact (left) and aproximated (right) solenoidal vector fields for Example 1, in a velocity mesh of 65 × 65 (a coarse mesh of 33 × 33).

Figure 5. Exact (left) and approximated (right) p for Example 1, in a coarse mesh of 33 × 33 (a velocity mesh of 65 × 65).

Table 2. Numerical results for Example 2. $\epsilon ={10}^{-4}$ .

method with and without preconditioning. We also solved a problem equivalent to (10) but with exact (not diagonal) stiffness matrix. A summary of the numerical results with the three algorithms is shown in Table 3, where we show the relative solenoidal error $e_r$ , average of divergence mdiv, as well as the number of iterations to get convergence, up to the given tolerance ( $\epsilon ={10}^{-12}$ , in the CG

Table 3. Numerical results for Example 3. $\epsilon ={10}^{-12}$ .

and PCG algorithms) for different mesh sizes. This table also contains the CPU time (sec) spent in each numerical experiment, and, in particular, in the last column the CPU time spent by solving the elliptic problem equivalent to (10). For this elliptic problem the results for $e_r$ and mdiv were not as good as those of the CG algorithms.

It is clear that the PCG-algorithm performs much better than the other algorithms. Another nice feature of the performance of the PCG-algorithm in this example is that the number of iterations is independent of the mesh size (7 iterations in each case). Also, the relative error between the computed and the exact solenoidal vector field is of the same order for the CG and PCG algorithms, with the largest difference occurring on the top boundary. For this example, we observe a loss of accuracy on the mean divergence when the PCG-algorithm is employed. Anyway, the mean divergence obtained with the PCG-algorithm is still very accurate, from a practical point of view, since we are using polynomial of degree 1 to approximate the unknowns.

5. Conclusions

We have applied an optimal preconditioner for the conjugate gradient algorithm, to solve the operator equation associated with the saddle-point formulation of the Helmholtz decomposition problem. Numerical results show that the new preconditioned conjugate gradient algorithm, and its stable approximation by a mixed finite element discretization, preserve the good properties of the non-preconditioned conjugate gradient algorithm, but it is much faster, and produces much better results than the methods based on the elliptic problem (10). Some additional nice properties of this algorithm are:

● It enforces very accurately mass conservation in the computed vector fields, considering that the approximation is second order.

● The number of iterations is reduced from several hundreds to less than 7 in the examples considered in this article. Also, the number of iterations of the PCG-algorithm is almost independent with respect to mesh refinement, while the number of iterations in the non-preconditioned algorithm doubles at each mesh refinement in most cases. Therefore, there is a substantial reduction in the computational time in both examples. However, this behavior must be corroborated in 3-D problems with adaptive meshes in complex domains, and with millions of degrees of freedom, in order to have a better idea of the performance of the method with realistic scenarios.

● It is not necessary to impose boundary conditions on the multiplier p, as it is done with the elliptic approaches. Furthermore, post-processing the scalar potential to find the vectorial field from the multiplier is not required, since the multiplier and the solenoidal vector field are found simultaneously by the algorithm.

The application of the methodology presented here to the more realistic three-dimensional case is an extension of the present work.

Acknowledgements

Sincere thanks to the members of JAMP for their professional attitude of high quality.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Agashe, S. (2021) New Formulas and Results for 3-Dimensional Vector Fields. Applied Mathematics, 12, 1058-1096. https://doi.org/10.4236/am.2021.1211069
[2]	Bhatia, H., Norgard, G., Pascucci, V. and Bremer, P. (2012) The Helmholtz-Hodge Decomposition—A Survey. IEEE Transactions on Visualization and Computer Graphics, 19, 1386-1404. https://doi.org/10.1109/TVCG.2012.316
[3]	Ahusborde, E., Azaez, M., Caltagirone, J.P., Gerritsma, M. and Lemoine, A. (2014) Discrete Hodge Helmholtz Decomposition. Monografas Matemáticas Garca de Galdeano, 39, 1-10.
[4]	Mansour, N.N., Kosovichev, A., Georgobiani, D., Wray, A. and Miesch, M. (2004) Turbulence Convection and Oscillations in the Sun. Proceedings of the SOHO14/GONG Workshop, New Haven, 12-16 July 2004, 164-171.
[5]	Akram, M. and Michel, V. (2010) Regularisation of the Helmholtz Decomposition and its Application to Geomagnetic Field Modelling. GEM—International Journal on Geomathematics, 1, 101-120. https://doi.org/10.1007/s13137-010-0001-y
[6]	Guo, Q., Mandal, M.K., Liu, G. and Kavanagh, K.M. (2006) Cardiac Video Analysis Using Hodge-Helmholtz Field Decomposition. Computers in Biology and Medicine, 36, 1-20. https://doi.org/10.1016/j.compbiomed.2004.06.011
[7]	Gao, H., Mandal, M., Guo, G. and Wan, J. (2010) Singular Point Detection Using Discrete Hodge Helmholtz Decomposition in Fingerprint Images. 2010 IEEE International Conference on Acoustics, Speech and Signal Processing, Dallas, 14-19 March 2010, 1094-1097. https://doi.org/10.1109/ICASSP.2010.5495348
[8]	Guo, Q., Mandal, M.K. and Li, M.Y. (2005) Efficient Hodge-Helmholtz Decomposition of Motion Fields. Pattern Recognition Letters, 26, 493-501. https://doi.org/10.1016/j.patrec.2004.08.008
[9]	Hinkle, J., Fletcher, P.T., Wang, B., Salter, B. and Joshi, S. (2009) 4D MAP Image Reconstruction Incorporating Organ Motion. In: Prince, J.L., Pham, D.L. and Myers, K.J., Eds., Information Processing in Medical Imaging, Springer, Heidelberg, 676-687. https://doi.org/10.1007/978-3-642-02498-6_56
[10]	Polthier, K. and Preu, E. (2003) Identifying Vector Fields Singularities Using a Discrete Hodge Decomposition. In: Hege, H.C. and Polthier, K., Eds., Visualization and Mathematics III, Springer, Heidelberg, 113-134. https://doi.org/10.1007/978-3-662-05105-4_6
[11]	Tong, Y., Lombeyda, S., Hirani, A. and Desbrun, M. (2003) Discrete Multiscale Vector Field Decomposition. ACM Transactions on Graphics, 22, 445-452. https://doi.org/10.1145/882262.882290
[12]	Foster, N. and Metaxas, D. (1997) Modeling the Motion of a Hot, Turbulent Gas. Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques, Los Angeles, 3-8 August 1997, 181-188. https://doi.org/10.1145/258734.258838
[13]	Ingber, M.S. and Kempka, S.N. (2001) A Galerkin Implementation of the Generalized Helmholtz Decomposition for Vorticity Formulations. Journal of Computational Physics, 169, 215-237. https://doi.org/10.1006/jcph.2001.6724
[14]	Deriaz, E. and Perrier, V. (2006) Divergence-Free and Curl-Free Wavelets in 2D and 3D: Application to Turbulent Flows. Journal of Turbulence, 7, 1-37. https://doi.org/10.1080/14685240500260547
[15]	Deriaz, E. and Perrier, V. (2008) Direct Numerical Simulation of Turbulence Using Divergence-Free Wavelets. Multiscale Modeling and Simulation, 7, 1101-1129. https://doi.org/10.1137/070701017
[16]	Deriaz, E. and Perrier, V. (2008) Orthogonal Helmholtz Decomposition in Arbitrary Dimension Using Divergence-Free and Curl-Free Wavelets. Applied and Computational Harmonic Analysis, 26, 249-269. https://doi.org/10.1016/j.acha.2008.06.001
[17]	Colin, F., Egli, R. and Lin, F. (2006) Computing a Null Divergence Velocity Field Using Smoothed Particle Hydrodynamics. Journal of Computational Physics, 217, 680-692. https://doi.org/10.1016/j.jcp.2006.01.021
[18]	Petronetto, F., Paiva, A., Lage, M., Tavares, G., Lopes, H. and Lewiner, T. (2010) Meshless Helmholtz-Hodge Decomposition. IEEE Transactions on Visualization and Computer Graphics, 16, 338-349. https://doi.org/10.1109/TVCG.2009.61
[19]	Stam, J. and Fiume, E. (1995) Depicting Fire and Other Gaseous Phenomena Using Diffusion Processes. Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Technique, Los Angeles, 6-11 August 1995, 129-136. https://doi.org/10.1145/218380.218430
[20]	Fuselier, E.J. (2006) Refined Error Estimates for Matrix-Valued Radial Basis Functions. Ph.D. Thesis, Texas A&M University, College Station.
[21]	Macedo, I. and Castro, R. (2010) Learning Divergence-Free and Curl-Free Vector Fields with Matrix-Valued Kernels. Technical Report, Instituto de Matemática Pura e Aplicada, Rio de Janeiro, Brazil.
[22]	Narcowich, F.J. and Ward, J.D. (1994) Generalized Hermite Interpolation via Matrix-Valued Conditionally Positive Definite Functions. Mathematics of Computation, 63, 661-687. https://doi.org/10.1090/S0025-5718-1994-1254147-6
[23]	López-López, J., Juárez, L.H. and Sandoval, M.L. (2016) Improving the Reconstruction of Vector Fields Using Mixed Finite Element Methods and Optimal Preconditioning. Numerical Methods for Partial Differential Equations, 32, 1137-1154. https://doi.org/10.1002/num.22044