1 Introduction↩︎

Meshless methods, i.e. numerical methods based on a set of points with no specific connectivity, have received great attention over the last two decades due to their ability to handle complex geometries, ability to adapt the points density in function of the local regularity and to support both Lagrangian and Eulerian formulations. Such an approach has been broken down into several methods such as Smooth Particles Hydrodynamics (SPH) [1], [2], Finite Pointset method (FPsM) [3], Finite Point method (FPtM) [4]–[6], or Generalized Finite Difference Method (GFDM) [7]–[9].

Finite Pointset and Finite Point methods use the entire neighborhood of points belonging to a compact ball, controlled by a threshold parameter included in the weighted function. Moreover, both rely on a moving least-squares method where the reconstruction is achieved at any point (not necessarily a point of the cloud). On the contrary, the Generalized Finite Difference Method depends on particular stencils (not all the points of the neighborhood are necessarily used), and a classical least-squares technique centered on the cloud points is used to derive relations between derivatives and the function values over the stencil. Consequently, the GFDM provides an additional advantage over FPsM and FPtM, the ability to play with the stencil choice and, by that way, obtain the best trade-off between stencil size and accuracy.

In [10], the authors take advantage of several properties inherent to the local representation in the GFDM. The choice of the stencil together with the parameterization of the weight function and the polynomial’s scaling parameter enable to optimize the quality of the representation of the derivatives while preserving small-sized stencils. The study was developed for the linear convection diffusion operator in two-dimensional space, where computational cost reduction and better accuracy for a given cloud of points are reported. The question then arises as to whether such gains are still present for non-linear problems. To this end, we conduct a new study on two very popular systems, the viscous Burgers’ equation [11] and the weakly compressible Navier-Stokes equations (Tait equation of state). We do not integrate any upwind techniques [12], [13] since we carry out the simulations with low Reynolds numbers.

2 Design of the numerical method↩︎

Let \(\Omega\) be a bounded open set of \(\mathbb{R}^2\) with boundary \(\partial \Omega\) and denote by \({\boldsymbol{x}}=(x_1,x_2)\) a generic point of \(\overline{\Omega}\). We denote by \(\phi\) any smooth function defined on \(\overline{\Omega}\) and sufficiently regular such that the first- and second-order derivatives are continuously bounded on \(\overline{\Omega}\). We restrict the presentation to the two-dimensional geometries for the sake of simplicity, but extension to higher dimensions is straightforward.

2.1 Cloud of points, approximations, and stencil↩︎

A cloud consists of a finite set of points \({\boldsymbol{x}}_i\in\overline{\Omega}\) characterized by the indices set \({\mathcal{G}}=\{1,\dots,I\}\) which we split into \({\mathring{\mathcal{G}}}\), the indices of the points that belong to the open set \(\Omega\) and \({\partial \mathcal{G}}\) the indices of points on the boundary \(\partial \Omega\). For any node \(i\in{\mathcal{G}}\), we denote by \({\mathcal{S}}_i\) the associated stencil constituted by a sublist of indices from \({\mathcal{G}}\) such that \(i\in{\mathcal{S}}_i\), while \(|{\mathcal{S}}_i|\) stands for the stencil size.

We introduce the multi-index notation \({\boldsymbol{\alpha}}=(\alpha_1,\alpha_2)\), with \(\alpha_1,\alpha_2\in \mathbb{N}_0\) with \(|{\boldsymbol{\alpha}}|=\alpha_1+\alpha_2\) and \({\boldsymbol{\alpha}}!=(\alpha_1!)(\alpha_2!)\), while \({\boldsymbol{x}}^{\boldsymbol{\alpha}}=(x_1)^{\alpha_1}(x_2)^{\alpha_2}\). For \(i\in{\mathcal{G}}\) we also define the monomial function of degree \({\boldsymbol{\alpha}}\) centred at point \({\boldsymbol{x}}_i\) with \[\pi_{\boldsymbol{\alpha}}({\boldsymbol{x}};{\boldsymbol{x}}_i,h_i)=\left ( \frac{{\boldsymbol{x}}-{\boldsymbol{x}}_i}{h_i}\right )^{\boldsymbol{\alpha}}\] where \(h_i>0\) is a scaling parameter.

Given \({\boldsymbol{\beta}}=(\beta_1,\beta_2)\), with \(\beta_1,\beta_2\in \mathbb{N}_0\), we associate the partial derivative \(\partial^{\boldsymbol{\beta}}\phi=\partial^{\beta_1}_{x_1}\,\partial^{\beta_2}_{x_2}\phi(x_1,x_2)\). Inequality \({\boldsymbol{\alpha}}\geq {\boldsymbol{\beta}}\) means that \(\alpha_1\geq \beta_1\) and \(\alpha_2\geq \beta_2\). In particular, if \({\boldsymbol{\alpha}}\geq {\boldsymbol{\beta}}\), one has \[\partial^{\boldsymbol{\beta}}\pi_{\boldsymbol{\alpha}}({\boldsymbol{x}};{\boldsymbol{x}}_i,h_i) =\frac{{\boldsymbol{\alpha}}!}{({\boldsymbol{\alpha}}-{\boldsymbol{\beta}})!}\,\frac{1}{(h_i)^{|{\boldsymbol{\beta}}|}}\,\pi_{{\boldsymbol{\alpha}}-{\boldsymbol{\beta}}}({\boldsymbol{x}};{\boldsymbol{x}}_i,h_i),\] while \(\partial^{\boldsymbol{\beta}}\pi_{\boldsymbol{\alpha}}({\boldsymbol{x}};{\boldsymbol{x}}_i,h_i)=0\), otherwise. Note that taking \({\boldsymbol{x}}={\boldsymbol{x}}_i\) provides the relations \[\label{eq::alpha95derivative95monomial} \partial^{\boldsymbol{\alpha}}\pi_{\boldsymbol{\alpha}}({\boldsymbol{x}}_i;{\boldsymbol{x}}_i,h_i)=\frac{{\boldsymbol{\alpha}}!}{(h_i)^{|{\boldsymbol{\alpha}}|}} \quad\textrm{ and }\quad \partial^{\boldsymbol{\beta}}\pi_\alpha({\boldsymbol{x}}_i;{\boldsymbol{x}}_i,h_i)=0, \textrm{ otherwise}.\tag{1}\]

2.2 Polynomial degree set and derivative order set↩︎

The polynomial degree set \(\mathcal{A}\). The space \(\mathbb{P}_i^{\mathcal{A}}\) is constituted by the span of the monomial family \[\mathbb{P}_i^{\mathcal{A}}=\textrm{span}\big\{\pi_{\boldsymbol{\alpha}}({\boldsymbol{x}};x_i,h_i),\;{\boldsymbol{\alpha}}\in {\mathcal{A}}\big \}\] we shall call the \(\mathcal{A}\)-polynomial space centred at node \(i\).

Consequently, any polynomial function of \(\mathbb{P}_i^{{\mathcal{A}}}\) is written under the form \[\pi({\boldsymbol{x}},{\boldsymbol{c}}_i;{\boldsymbol{x}}_i)=\pi({\boldsymbol{x}},{\boldsymbol{c}}_i;{\boldsymbol{x}}_i,h_i,{\mathcal{A}})=\sum_{{\boldsymbol{\alpha}}\in{\mathcal{A}}} c_{i,{\boldsymbol{\alpha}}}\pi_{\boldsymbol{\alpha}}({\boldsymbol{x}};{\boldsymbol{x}}_i,h_i),\] where the vector \({\boldsymbol{c}}_i\) collects the entries \(\big (c_{i,{\boldsymbol{\alpha}}}\big)_{{\boldsymbol{\alpha}}\in{\mathcal{A}}}\). The set \({\mathcal{B}}\) contains the list of derivatives \({\boldsymbol{\beta}}\) we aim at discretizing via differential operators. For example, \(\big \{\partial^0_{x_1}\partial^0_{x_2},\, \partial^1_{x_1}\partial^0_{x_2},\, \partial^0_{x_1}\partial^1_{x_2},\,\partial^1_{x_1}\partial^1_{x_2}\big \}\iff {\mathcal{B}}=\big \{(0,0),(1,0),(0,1),(1,1)\big \}\).

Finally, for any regular function \(\phi\), we denote \(\phi_i\approx \phi({\boldsymbol{x}}_i)\) an approximation of \(\phi\) at point \({\boldsymbol{x}}_i\) and for any partial derivative \(\partial^{\boldsymbol{\beta}}\), we use the notation \(\phi_i^{\boldsymbol{\beta}}\approx \partial^{\boldsymbol{\beta}}\phi({\boldsymbol{x}}_i)\). We then organize the data following two vectors:

• \(\Phi_{{\mathcal{S}}_i}=(\phi_j)_{j\in{\mathcal{S}}_i}\) the approximation vector associated with the stencil \({\mathcal{S}}_i\);

• \(\Phi=\Phi_{{\mathcal{G}}}=(\phi_i)_{i\in{\mathcal{G}}}\) the vector that gathers all the approximations associated to the cloud \({\mathcal{G}}\).

2.3 Interpolation and differential operators’ discretization↩︎

Given a stencil \({\mathcal{S}}_i\) associated to a node \(i\in{\mathcal{G}}\), \(\Phi_{{\mathcal{S}}_i}\) the vector of approximations over the stencil, and \(\mathbb{P}^{\mathcal{A}}_i\) the polynomial subspace, we define the following functional for any \(\pi\in\mathbb{P}_i^{\mathcal{A}}\), \[E({\boldsymbol{c}}_i;{\mathcal{S}}_i,\Phi_{{\mathcal{S}}_i},{\boldsymbol{w}}_i,{\mathcal{A}},h_i)=\sum_{j\in{\mathcal{S}}_i} w_{ij}\Big [\pi({\boldsymbol{x}}_j,{\boldsymbol{c}}_i;{\boldsymbol{x}}_i,h_i,{\mathcal{A}})-\phi_j\Big ]^2,\] where \({\boldsymbol{w}}_i=(w_{ij})_{j\in{\mathcal{S}}_i}\) contains the weights associated to each point of the stencil \({\mathcal{S}}_i\). A usual choice is \(w_{ij}=\rho(\Vert {\boldsymbol{x}}_i-{\boldsymbol{x}}_j\Vert;\tau_i)\), where \(\rho(s;\tau)\) is a decreasing function of \(s\in [0,+\infty[\) controlled by a scaling parameter \(\tau\). For instance, the Gaussian weights read \[w_{ij}=\rho(\Vert {\boldsymbol{x}}_j-{\boldsymbol{x}}_i\Vert;\tau_i)=\exp \left( - \frac{\Vert {\boldsymbol{x}}_j-{\boldsymbol{x}}_i\Vert^2 }{\tau^2_i} \right).\]

We define the \(|{\mathcal{A}}|\times|{\mathcal{S}}_i|\) matrix \(\displaystyle M_i=\big (m_{i;{\boldsymbol{\alpha}},j}\big )^{{\boldsymbol{\alpha}}\in{\mathcal{A}}}_{j\in{\mathcal{S}}_i}\) with \[m_{i;{\boldsymbol{\alpha}},j}=\pi_{\boldsymbol{\alpha}}({\boldsymbol{x}}_j;{\boldsymbol{x}}_i,h_i)= \left ( \frac{{\boldsymbol{x}}_j-{\boldsymbol{x}}_i}{h_i}\right )^{\boldsymbol{\alpha}},\quad {\boldsymbol{\alpha}}\in{\mathcal{A}},\, j\in{\mathcal{S}}_i.\] Using the Euclidean norm, the functional rewrites under the form \[E_i({\boldsymbol{c}}_i;{\mathcal{S}}_i,\Phi_{{\mathcal{S}}_i},{\boldsymbol{w}}_i,{\mathcal{A}},h_i)= \Vert W_i(M_i{\boldsymbol{c}}_i-\Phi_{{\mathcal{S}}_i})\Vert^2,\] with \(W_i=\textrm{diag}({\boldsymbol{w}}_{i})\) the \(|{\mathcal{S}}_i|\times |{\mathcal{S}}_i|\) diagonal matrix of the weights.

Minimization of the functional with respect to vector \({\boldsymbol{c}}_i\) leads to the linear system \[(M_i)^ {{\textrm t}} W_i M_i\, {\boldsymbol{c}}_i=(M_i)^ {{\textrm t}} W_i\Phi_{{\mathcal{S}}_i}.\] Assuming that the choice of \({\mathcal{A}}\) and \({\mathcal{S}}_i\) provides a matrix \(M_i\) that satisfies the maximal rank property, then \(M_i^{\dagger}=(M_i)^ {{\textrm t}} W_i M_i\) is a non-singular \(|{\mathcal{A}}|\times|{\mathcal{A}}|\) matrix and the polynomial coefficients are given by \({\boldsymbol{c}}_i=M_i^\ddag \Phi_{{\mathcal{S}}_i}\) with \[M_i^\ddag =\big (M_i^\dagger\big )^{-1}(M_i)^ {{\textrm t}} W_i.\] Consequently, we deduce the coefficient of a \(\mathcal{A}\)-polynomial approximation given by \[\widetilde{\phi}_i({\boldsymbol{x}})=\phi_i({\boldsymbol{x}},{\boldsymbol{c}}_i;{\boldsymbol{x}}_i,{\mathcal{S}}_i,h_i,\mathcal{A})=\sum_{\alpha\in \mathcal{A}}c_{i,{\boldsymbol{\alpha}}} \pi_{\boldsymbol{\alpha}}({\boldsymbol{x}};{\boldsymbol{x}}_i,h_i).\] Under the condition \(\mathcal{B}\subset \mathcal{A}\), we get \[\partial^{\boldsymbol{\beta}}\pi_\alpha({\boldsymbol{x}}_i)=\delta_{{\boldsymbol{\alpha}},{\boldsymbol{\beta}}}\frac{{\boldsymbol{\alpha}}!}{(h_i)^{\boldsymbol{\alpha}}}.\] and deduce \[\label{approximation95discretization} \partial^\beta \widetilde{\phi}_i({\boldsymbol{x}}_i)=\frac{{\boldsymbol{\beta}}!}{(h_i)^{\boldsymbol{\beta}}}c_{i,{\boldsymbol{\beta}}} =\frac{{\boldsymbol{\beta}}!}{(h_i)^{|{\boldsymbol{\beta}}|}}\, \sum_{j\in{\mathcal{S}}_i} m^\ddag_{i;{\boldsymbol{\beta}},j} \phi_j =\sum_{j\in{\mathcal{S}}_i} a_{i;{\boldsymbol{\beta}},j} \phi_j, \textrm{ with } a_{i;{\boldsymbol{\beta}},j}=\frac{{\boldsymbol{\beta}}!}{(h_i)^{|{\boldsymbol{\beta}}|}}\,m^\ddag_{i;{\boldsymbol{\beta}},j}.\tag{2}\] Relation 2 provides the discretization of the partial derivatives, in terms of the stencil-node values, which we use in the physical equations.

2.4 Discretization optimization↩︎

Given a point \({\boldsymbol{x}}_i\in{\mathcal{G}}\), matrix \(M_i^\ddag=M_i^\ddag(h_i,\tau_i,{\mathcal{S}}_i)\) depends on the polynomial scaling factor \(h_i\), the weights scaling factor \(\tau_i\), and the stencil \({\mathcal{S}}_i\). In [10] we develop an optimization procedure based on the condition number of matrix \(M_i^\dagger\) since it corresponds to the non-singular matrix we have to inverse to determine the derivatives. We then introduce a quantification of the polynomial reconstruction efficiency by computing the condition number \(\chi_2(M^\dagger_i)\) and we aim at determining the parameters that minimize it. In other words, we seek \[(h_i^\star,\tau_i^\star,{\mathcal{S}}^\star_i)=\arg \min \, \chi_2(M_i^\dagger(h_i,\tau_i,{\mathcal{S}}_i))\] under the constraint \(|{\mathcal{S}}_i|\geq |{\boldsymbol{\alpha}}|\), since we need more points (i.e., equations) than unknown polynomial coefficients.

Optimization over a discrete set \({\mathcal{S}}\) and continuous sets \(]0,+\infty[\times ]0,+\infty[\) is a difficult task, so the authors combine a heuristic approach for the stencil choice along with a gradient method for the continuous parameters. In this way, we obtain a rather small stencil (\(|{\mathcal{A}}|\leq |{\mathcal{S}}_i|\leq 1.5 |{\mathcal{A}}|\)) with a dramatic cut in the condition number in comparison with an unsophisticated stencil consisting of the closest points to the reference point \({\boldsymbol{x}}_i\) (about twice the size of \({\mathcal{A}}\) to produce a non-singular matrix \(M_i^\dagger\)).

Since we report excellent orders of accuracy for linear PDE problems on complex domains (curved boundary domains) in [10], we aim to extend the same framework to non-linear PDE problems and systems by benchmarking several important non-linear problems and assessing its efficiency in terms of both accuracy and stability.

3 The two-dimensional Burgers’ equations↩︎

The main goal is to assess the accuracy improvement achieved by the optimized stencil, as well as the robustness of the method. The linear case has been studied in [10] where we showed the advantages and gains of the optimization, but the question arises when we are dealing with non-linear problems.

3.1 Modeling and linearization↩︎

We consider the two-dimensional stationary forced Burgers’ equations with viscous term, \[\begin{align} u_1\,\partial_{x_1}u_1+u_2\,\partial_{x_2}u_1-\varepsilon \big( \partial_{x_1}^{2}u_1+\partial_{x_2}^{2}u_1\big) & =f_{1},\quad \textrm{ in } \Omega , \\ u_1\,\partial_{x_1}u_2+u_2\,\partial_{x_2}u_2-\varepsilon \big( \partial_{x_1}^{2}u_2+\partial_{x_2}^{2}u_2\big) & =f_{2},\quad \textrm{ in } \Omega , \end{align}\] where \({\boldsymbol{u}}=(u_1,u_2)\) is the unknown velocity, \({\boldsymbol{f}}=(f_1,f_2)\) the source term, and \(\varepsilon > 0\) the viscosity parameter. The general Robin boundary conditions are given by \[\mu_1({\boldsymbol{x}}_i)\, u_1+\nu_1({\boldsymbol{x}}_i) \nabla u_1\cdot {\boldsymbol{n}}=g_1({\boldsymbol{x}}_i),\quad \mu_2({\boldsymbol{x}}_i)\, u_2+\nu_2({\boldsymbol{x}}_i) \nabla u_2\cdot {\boldsymbol{n}}=g_2({\boldsymbol{x}}_i), \quad \textrm{ on } \partial \Omega,\] where \({\boldsymbol{\mu}}=(\mu_1,\mu_2)\), \({\boldsymbol{\nu}}=(\nu_1,\nu_2)\), and \({\boldsymbol{g}}=(g_1,g_2)\) are given functions and \({\boldsymbol{n}}={\boldsymbol{n}}({\boldsymbol{x}})\) is the outward unit normal vector. Notice that the velocity components are not coupled on the boundary. In this context, we define for any node \(i\in{\mathcal{G}}\) the label \(\ell(i)\) that indicates the physical equation associated with the node. For example, \(\ell(i)=1\) means that one has to satisfy the inner equations, while \(\ell(i)=2\) corresponds to the case where the boundary equations hold.

To solve the non-linear system, we use a simple fixed-point method based on a local linearization where we build a sequence \({\boldsymbol{u}}^{[k]}=(u_1^{[k]},u_2^{[k]})\), \(k=0,1,\dots\), of approximations of the solution \({\boldsymbol{u}}\) on \(\overline{\Omega}\). To this end, assuming \((u_1^{[k]},u_2^{[k]})\) is known, we seek a solution \((u_1^{[k+1]},u_2^{[k+1]})\) of the linearized problem \[\begin{align} u_1^{[k]}\,\partial_{x_1}u_1^{[k+1]}+u_2^{[k]}\,\partial_{x_2}u_1^{[k+1]}-\varepsilon\Big (\partial_{x_1}^{2}u_1^{[k+1]}+\partial_{x_2}^{2}u_1^{[k+1]}\Big )&=f_{1}, \\ u_1^{[k]}\,\partial_{x_1}u_2^{[k+1]}+u_2^{[k]}\,\partial_{x_2}u_2^{[k+1]}-\varepsilon\Big (\partial_{x_1}^{2}u_2^{[k+1]}+\partial_{x_2}^{2}u_2^{[k+1]}\Big )&=f_{2}. \end{align}\] Note that the problem is decoupled and one solves for the velocities \(u_1\) and \(u_2\) independently. The iterative procedure stops when both \(\Vert u_1^{[k+1]}-u_1^{[k]}\Vert\) and \(\Vert u_2^{[k+1]}-u_2^{[k]}\Vert\) are lower than a given tolerance. Convergence of the fixed-point is not guaranteed, in particular for large Péclet numbers, since the stencil choice does not take into account the upwind direction, leading to a centered scheme for the convective term. Therefore, benchmarks have been carried out with a relatively low Péclet number to focus on the accuracy and convergence order of the numerical scheme.

3.2 Discretization and assembly↩︎

We denote by \(U_1\) and \(U_2\) the vectors of \(\mathbb{R}^I\) with the numerical solution of the non-linear system, that is, \[U_1=(u_{1,i})_{i\in{\mathcal{G}}}, \quad U_2=(u_{2,i})_{i\in{\mathcal{G}}}.\] Moreover, for any point \({\boldsymbol{x}}_i\), \(i\in{\mathcal{G}}\), and \(\phi=u_1,u_2\), we have \(\phi({\boldsymbol{x}}_i)\approx\phi_i^{(0,0)}=\phi_i\), and we use relations 2 for the discrete derivatives, that is \[\begin{align} &\displaystyle \partial_{x_1}\phi({\boldsymbol{x}}_i)\approx\phi_i^{(1,0)}= \sum_{j\in{\mathcal{S}}_i} a_{i;(1,0),j} \phi_j, &\displaystyle \partial_{x_2}\phi({\boldsymbol{x}}_i)\approx\phi_i^{(0,1)}= \sum_{j\in{\mathcal{S}}_i} a_{i;(0,1),j} \phi_j,\tag{3}\\ &\displaystyle \partial^2_{x_1}\phi({\boldsymbol{x}}_i)\approx\phi_i^{(2,0)}= \sum_{j\in{\mathcal{S}}_i} a_{i;(2,0),j} \phi_j, &\displaystyle \partial^2_{x_2}\phi({\boldsymbol{x}}_i)\approx\phi_i^{(0,2)}= \sum_{j\in{\mathcal{S}}_i} a_{i;(0,2),j} \phi_j.\tag{4} \end{align}\] Since the simulation is carried out using a simple fixed-point method, we construct a sequence of approximations using the following iterative procedure: given (\(U_1^{[k]},U_2^{[k]}\)), we compute \((U_1^{[k+1]},U_2^{[k+1]})\), the approximation of the associated continuous linearized problem at iteration \(k+1\): \[\begin{align} u_{1,i}^{[k],(0,0)}\,u_{1,i}^{[k+1],(1,0)}+u_{2,i}^{[k],(0,0)}\,u_{1,i}^{[k+1],(0,1)}-\varepsilon \left( u_{1,i}^{[k+1],(2,0)}+u_{1,i}^{[k+1],(0,2)}\right ) & =f_{1}({\boldsymbol{x}}_i), \\ u_{1,i}^{[k],(0,0)}\,u_{2,i}^{[k+1],(1,0)}+u_{2,i}^{[k],(0,0)}\,u_{2,i}^{[k+1],(0,1)}-\varepsilon \left(u_{2,i}^{[k+1],(2,0)}+u_{2,i}^{[k+1],(0,2)}\right ) & =f_{2}({\boldsymbol{x}}_i). \end{align}\] Replacing the derivatives by the local representation given by 3 and 4 , we then assemble two independent linear systems \(A^{[k]}_1U_1=b_1\) and \(A^{[k]}_2U_2=b_2\). For example, the assembled matrix corresponding to vector \(U_1\) is obtained as follows:

• if \(\ell(i)=1\), then for \(j\in{\mathcal{S}}_i\) we set \[A^{[k]}_1[i,j]=u_{1,j}^{[k]}\,a_{i;(1,0),j}+u_{2,j}^{[k]}\,a_{i;(0,1),j}-\varepsilon\Big(a_{i;(2,0),j}+a_{i;(0,2),j} \Big )\] and \(A^{[k]}_1[i,j]=0\) otherwise. The right-hand side is simply \(b_1[i]=f_1({\boldsymbol{x}}_i)\).

• if \(\ell(i)=2\), then for \(j\in{\mathcal{S}}_i\) the discrete Robin condition provides the following entries \[A^{[k]}_1[i,j]=\delta_{i,j}\mu_1({\boldsymbol{x}}_i)+\nu_1({\boldsymbol{x}}_i) (a_{i;(1,0),j}n_{i;1}+a_{i;(0,1),j}n_{i;2}),\] and \(A^{[k]}_1[i,j]=0\) otherwise. The right-hand side is given by \(b_1[i]=g_1({\boldsymbol{x}}_i)\). Note that whenever the Dirichlet conditions hold, one has \(A^{[k]}_1[i,j]=\delta_{i,j}\) and \(b_1[i]=u_{D_1}({\boldsymbol{x}}_i)\).

We update the numerical approximation by computing \(U_1^{[k+1]}\) and \(U_2^{[k+1]}\), solutions of \(A^{[k]}_1\,U_1^{[k+1]}=b_1\) and \(A^{[k]}_2\,U_2^{[k+1]}=b_2\), respectively.

3.3 Test case 3.1↩︎

In this first test case devoted to the two-dimensional Burgers’ equations, we prescribe the Dirichlet boundary conditions for both horizontal and vertical velocities, that is \({\boldsymbol{\mu}}({\boldsymbol{x}}_i)=\mathbf{1}\), \({\boldsymbol{\nu}}({\boldsymbol{x}}_i)=\mathbf{0}\), for all \(i\in{\partial \mathcal{G}}\). We further take \(\varepsilon=10^{-2}\) and consider the following irregular domain, \[\Omega =\left\{ {\boldsymbol{x}}\in \left] -1,1\right[ \times \left] -1,1\right[ :x_1^{2}+x_2^{2}>0.3^{2}\wedge |x_1x_2|<0.3\right\}.\] In order to perform a convergence analysis we set \({\boldsymbol{f}}=(f_1,f_2)\) and \({\boldsymbol{u}}_D=(u_{D_1},u_{D_2})\) such that \[u_1({\boldsymbol{x}})=-\frac{\varepsilon x_2}{x_1^{2}+x_2^{2}},\qquad u_2({\boldsymbol{x}})=\frac{\varepsilon x_1}{x_1^{2}+x_2^{2}}\] is the unique solution of the boundary value problem (BVP). We note that this solution is invariant to rotation with respect to the origin of the coordinate system. The resulting expression for \((f_1,f_2)\) is \[f_{1}({\boldsymbol{x}})=-\frac{\varepsilon ^{2}x_1}{\left( x_1^{2}+x_2^{2}\right) ^{2}},\qquad f_{2}({\boldsymbol{x}})= \frac{\varepsilon ^{2}x_2}{\left( x_1^{2}+x_2^{2}\right) ^{2}}.\] The initial velocity for the fixed-point method, (\(U_1^{[0]},U_2^{[0]}\)), is taken equal to the average value over the nodes located on \(\partial \Omega\). The geometrical setting and the exact solution of the BVP are illustrated in Figure 1.

Numerical simulations are carried out with \(622\), \(1228\), \(2458\), and \(4886\) nodes, obtained with a Delaunay triangulation procedure using the Gmsh\(\textsuperscript{\copyright}\) package. We note that the number of nodes on the inner and outer boundaries is almost the same. Figure 2 illustrates the different clouds used in the present test case.

We present the numerical solution error (using \(L^{\infty}\)- and \(L^{1}\)-norms) and the corresponding average convergence order (ACO) for \(p=2,\ldots,6\) in Figure 3.

Looking at ACO, we conclude that the optimal convergence order is obtained (or even exceeded) for both the \(L^{\infty}\)- and the \(L^1\)-errors, independently of the approximation degree. We also report that the ratio \(L^{\infty}\)-error / \(L^1\)-error obtained for \(p=2,\ldots,6\) (\(4.9 \pm 0.5\), \(4.1 \pm 0.2\), \(4.2 \pm 0.4\), \(5.9 \pm 0.9\), and \(9.8 \pm 2.0\), respectively) is well bounded, indicating the good quality of the numerical solution across the computational domain for all values of \(p\).

3.4 Test case 3.2↩︎

We consider once again the two-dimensional Burgers’ equations supplemented by Dirichlet boundary conditions over a non-polygonal domain expressed in polar coordinates with respect to the origin of the coordinate system (see Figure 4): \[\Omega =\left\{ (r,\theta)\in \mathbb{R}^{+} \times \left[ 0,2\pi\right[ : 0.3 \leq r \leq 0.95+0.05\sin(12\theta+\pi/2)\right\}.\]

In order to perform a convergence analysis, we set \({\boldsymbol{f}}\) and \({\boldsymbol{u}}_D\) such that the unique solution of the BVP is given by \[u_1(r,\theta)=-q(r,\varepsilon)\sin\theta, \qquad u_2(r,\theta)=q(r,\varepsilon)\cos\theta,\] with \[q(r,\varepsilon )=\frac{r_4}{\varepsilon }\exp \left[ \frac{r_{3}\varepsilon }{\left( r-r_{1}\right) \left( r-r_{2}\right) }\right],\] where we adopt \(\varepsilon=10^{-2}\), while \(r_1\), \(r_2\), \(r_3\), and \(r_4\) are free parameters we shall set in the sequel. The solution of the BVP is invariant to rotation with respect to the origin of the coordinate system and we obtain \[\begin{align} f_{1}( r,\theta) =-\frac{q^{2}}{r}\cos \theta +\varepsilon \left( \frac{\partial ^{2}q}{\partial r^{2}}+\frac{1}{r\mathstrut }\frac{\partial q}{\partial r}-\frac{q}{r^{2}}\right) \sin \theta, \\ f_{2}(r,\theta) =-\frac{q^{2}}{r}\sin \theta -\varepsilon \left( \frac{\partial ^{2}q}{\partial r^{2}}+\frac{1}{r\mathstrut }\frac{\partial q}{\partial r}-\frac{q}{r^{2}}\right) \cos \theta. \end{align}\]

3.4.0.1 Test case 3.2a

We start the numerical simulations by setting \(r_1=3\), \(r_2=0.225\), \(r_3=100\), and \(r_4=2\times 10^{-4}\). The corresponding shape for the function \(q\) (leading to a shallow slope curve) and the exact solution of the BVP are depicted in Figure 5.

For the numerical simulations, we use four grids with \(674\), \(1316\), \(2533\), and \(5118\) nodes, obtained with a Delaunay triangulation procedure. The number of nodes in the inner boundary doubles that of its outer counterpart. Figure 6 illustrates the different grids used in the test case 3.2a numerical simulations. The nonlinearity is addressed using a fixed-point approach, as described previously.

In Figure 7, we present the numerical solution error (using \(L^{\infty}\)- and \(L^{1}\)-norms) and the corresponding ACO for \(p=2,\ldots,6\).

From the analysis of the results obtained for the ACO, we conclude that the optimal convergence order is obtained (or even exceeded) for both the \(L^{\infty}\)- and the \(L^1\)-errors, as long as the approximation degree is even. When \(p\) is odd, the ACO obtained equals the one obtained for the \(p-1\) case, a feature already reported in our previous work.

We also report that the ratio \(L^{\infty}\)-error / \(L^1\)-error obtained for \(p=2,\ldots,6\) (\(5.1 \pm 0.8\), \(3.2 \pm 0.2\), \(5.3 \pm 0.7\), \(4.0 \pm 0.7\), and \(9.9 \pm 4.0\), respectively) remains moderate and very similar to the values obtained in case 3.1, indicating the good quality of the numerical solution for the whole computational domain for all values of \(p\).

3.4.0.2 Test case 3.2b

We now consider a solution that has an increased gradient near the inner boundary by setting \(r_1=4\), \(r_2=0.2775\), \(r_3=50\), and \(r_4=2.8\times 10^{-4}\). The corresponding shape of function \(q\) and the exact solution of the BVP are illustrated in Figure 8.

For the numerical simulations, we use four grids with \(667\), \(1317\), \(2521\), and \(5178\) nodes, obtained with a Delaunay triangulation procedure. In order to account for the increased gradient near the inner boundary, the number of nodes in the inner boundary is now four times the number of nodes in its outer counterpart. We note that the total number of nodes in the four grids is very similar to the one in test case 3.2a. Figure 9 illustrates the different grids used in the test case 3.2b simulations.

The results obtained for the numerical solution error and convergence order for \(p=2,\ldots,6\) are presented in Figure 10.

We conclude that even though the errors are higher than those obtained in test case 3.2a, increasing the gradient of function \(q\) near the inner boundary does not impact the convergence order, since it follows the same trends already described for test case 3.2a. Likewise, the results observed for the ratio \(L^{\infty}\)-error / \(L^1\)-error for \(p=2,\ldots,6\) (\(5.1 \pm 0.8\), \(4.3 \pm 0.5\), \(4.7 \pm 1.5\), \(4.9 \pm 1.1\), \(7.0 \pm 1.2\), respectively) are also in line with the good values obtained for the previous test cases.

3.4.0.3 Test case 3.2c

Finally, setting \(r_1=2\), \(r_2=0.295\), \(r_3=5\), and \(r_4=3.45\times 10^{-4}\), we consider a solution that is even steeper near the inner boundary, mimicking a boundary layer scenario. The corresponding shape for function \(q\) and the exact solution of the BVP are shown in Figure 11.

In order to account for the rather steep gradient near the inner boundary, the number of nodes in the inner boundary is set to ten times the number of nodes in its outer counterpart. The density of nodes in the outer region has to be kept similar to the one used in test case 3.2b in order to ensure the convergence of the fixed point method for all values of \(p\). Therefore, the number of nodes has to be increased with respect to the test cases 3.2a and 3.2b. Hence, for the numerical simulations, we use four grids with \(2711\), \(3723\), \(5159\), and \(7169\) nodes, obtained with a Delaunay triangulation procedure. Figure 12 illustrates the four grids used in the test case 3.2c simulations.

The results obtained for the numerical solution error and the ACO for \(p=2,\ldots,6\) are presented in Figure 13.

The results show that the use of a more demanding function \(q\) inevitably leads to higher errors, but the convergence order remains essentially the same as in test cases 3.2a and 3.2b, reflecting the good ability of the numerical scheme to handle a boundary layer-like scenario. On the other hand, the results observed for the ratio \(L^{\infty}\)-error / \(L^1\)-error for \(p=2,\ldots,6\) (\(4.2 \pm 0.6\), \(4.1 \pm 0.4\), \(4.8 \pm 1.6\), \(3.2 \pm 0.5\), \(10.6 \pm 3.2\), respectively) support the good quality of the numerical solution across the computational domain for all values of \(p\).

3.5 Test case 3.3↩︎

For this last test case concerning the two-dimensional Burgers’ equations, we aim to test the ability of the numerical scheme to address Neumann boundary conditions in the presence of complex geometries. As in the previous test case, we take \(\varepsilon=10^{-2}\). We now consider a 2D nozzle-like domain given by

\[\Omega =\left\{ {\boldsymbol{x}}\in \mathbb{R}^{2}:0.5<x_1<1,\;-h(x_1)<x_2<h(x_1)\right\} ,\] where \[\begin{align} \Gamma_{t}& =\left\{ {\boldsymbol{x}}\in \mathbb{R}^{2}:0.5<x_1<1,\;x_2=h(x_1)\right\} , \\ \Gamma_{b}& =\left\{ {\boldsymbol{x}}\in \mathbb{R}^{2}:0.5<x_1<1,\;x_2=-h(x_1)\right\} , \\ \Gamma_{l}& =\left\{ {\boldsymbol{x}}\in \mathbb{R}^{2}:x_1=0.5,\;-h(x_1)\leq x_2\leq h(x_1)\right\} , \\ \Gamma_{r}& =\left\{ {\boldsymbol{x}}\in \mathbb{R}^{2}:x_1=1,\;-h(x_1)\leq x_2\leq h(x_1)\right\} , \end{align}\]

with \(h(x_1)=0.025\arccos(0.1\cos40x_1)\). The solution of the BVP satisfies homogeneous Neumann conditions on the upper and lower boundaries, \(\Gamma_{t}\) and \(\Gamma_{b}\), and Dirichlet conditions on the left and right boundaries, \(\Gamma_{l}\) and \(\Gamma_{r}\), that is, \[\begin{array}{l} {\boldsymbol{u}}({\boldsymbol{x}})={\boldsymbol{u}}_{D}({\boldsymbol{x}}),\text{\quad on } \Gamma_{l}\cup \Gamma _{r},\smallskip \\ \partial _{\mathbf{n}}{{\boldsymbol{u}}}({\boldsymbol{x}}) = \mathbf{0},\text{\quad on } \Gamma _{t}\cup \Gamma _{b}, \end{array}\] where \(\mathbf{n}\) is the outer unit normal vector. Therefore, we set \({\boldsymbol{\mu}}({\boldsymbol{x}}_i)=\mathbf{1}\) and \({\boldsymbol{\nu}}({\boldsymbol{x}}_i)=\mathbf{0}\), for all \({\boldsymbol{x}}_i \in\Gamma_{l} \cup \Gamma_{r}\), and \({\boldsymbol{\mu}}({\boldsymbol{x}}_i)=\mathbf{0}\) and \({\boldsymbol{\nu}}({\boldsymbol{x}}_i)=\mathbf{1}\), for all \({\boldsymbol{x}}_i \in\Gamma_{t} \cup \Gamma_{b}\). In order to perform a convergence analysis we set \({\boldsymbol{f}}\) and \({\boldsymbol{u}}_D\) such that \[u_1({\boldsymbol{x}})=\frac{1}{5}(1+\cos40x_1\cos40x_2),\text{\quad }u_2({\boldsymbol{x}})=\frac{1}{5}(1-\cos40x_1\cos40x_2),\] is the unique solution of the BVP. The choice of function \(h\) ensures homogeneous Neumann conditions hold on \(\Gamma _{t}\) and \(\Gamma_{b}\) for this particular choice of the solution of the BVP (cf. [14]).

The geometrical setting and the exact solution of the BVP are illustrated in Figure 14.

For the numerical simulations, we use four grids with \(625\), \(1264\), \(2505\), and \(5008\) nodes, obtained using a Delaunay triangulation procedure. Figure 15 illustrates the different grids used in test case 3.3 numerical simulations.

The results obtained for the numerical solution error and convergence order for \(p=2,\ldots,6\) are presented in Figure 16.

From the analysis of the results, we conclude that optimal convergence orders are obtained for the two norms used to evaluate the numerical solution error when the polynomial approximation degree is even. As in test case 3.2, the results obtained for \(p=3\) and \(p=5\), in particular the ACO, are rather similar to those found for \(p=2\) and \(p=4\), respectively. The good quality of the numerical solution across the computational domain is also confirmed by the low values obtained for the ratio \(L^{\infty}\)-error / \(L^1\)-error (\(3.9 \pm 0.8\), \(2.6 \pm 0.3\), \(3.4 \pm 0.9\), \(3.1 \pm 1.1\), \(3.7 \pm 0.5\), for \(p=2,\ldots,6\), respectively), which remains below \(4\).

4 Conclusion↩︎

The general finite difference method, by its nature, enables one to wisely choose the nodes of the stencil attached to a generic node \(i\) while other popular meshless methods, such as FPsM or FPtM, always consider the nodes belonging to the compact ball given by a user-dependent radius. Such a degree of freedom provides both accuracy and stability while maintaining a small stencil size. As a continuation of our work on the scalar linear PDE, we prove that the strategy works very well for the non-linear steady-state problem using a fixed-point technique based on a local linearization. Only low local Péclet/Reynolds-number discrete problems were considered to assess the optimization algorithm initially based on an isotropic stencil. The next stage will be to introduce some upwind bias, both for the stencil and the weights, to tackle convection-dominated problems and even hyperbolic systems with smooth solutions.

Acknowledgments↩︎

The authors Stéphane Clain and Jorge Figueiredo were funded by national funds from FCT - Fundação para a Ciência e a Tecnologia, I.P., within the context of the project HiPMos4C setting with reference 2023.16854.ICDT2023.16854.ICDT, https://doi.org/10.54499/2023.16854.ICDT.
Stéphane Clain would like to acknowledge the support of Centro de Matemática of the Coimbra University (CMUC), funded by the Portuguese Foundation for Science and Technology, under the references UIDB/00324/2020 (https://doi.org/10.54499/UIDB/00324/2020) and UIDP/00324/2020 (https://doi.org/10.54499/UIDP/00324/2020).
Jorge Figueiredo was supported by national funds through FCT - Fundação para a Ciência e a Tecnologia, under the Program Contract of the Centre of Mathematics of the University of Minho (CMAT/UM), UID/00013/2025 (https://doi.org/10.54499/UID/00013/2025).

Conflicts of Interest↩︎

The authors declare that there are no conflicts of interest.

References↩︎