In physical modelling, systems of partial differential equations (PDEs) are used to describe the dynamical properties of many natural phenomena. Moreover, the solution of such systems is often of interest to engineers and scientists. However, due to their complexity, they almost never have analytical solutions, and need to be treated numerically, leading to a numerical solution. In general, PDE problems are often solved using one of the following three methods: the finite volume method (FVM), the finite element method (FEM) and the finite difference method (FDM). Recently, however, a generalised formulation of FDM, the radial basis function-generated finite differences (RBF-FD)[tolstykh2000using, tolstykh2003using], has become increasingly popular. This is mainly because RBF-FD is a variant of the mesh-free methods [belytschko1996meshless], i.e. the method can operate on scattered nodes, unlike the previously mentioned mesh-based methods.
In the context of RBF-FD, linear differential operators are approximated over a set of RBFs augmented with monomials. Augmentation is necessary to ensure convergent and stable behaviour of the method [bayona2017role, flyer2016]. Additionally, it also enables a direct control over the order of the approximation method, as it corresponds to the highest order used in the approximation basis.
Nevertheless, after the numerical solution is obtained, scientists are often confronted with the difficulty of validating it. For that reason, researchers proposed error indicators [9597066, carstensen2005explicit] to identify problematic areas with a high error of the numerical solution. In practise, different adaptive numerical methods are then applied to these areas [SEGETH20101589] ensuring a finer local field description (h-adaptivity) or higher polynomial degree approximations (p-adaptivity), effectively improving the accuracy of numerical solution.
In this paper, we present an a posteriori error indicator that measures the error of an implicit solution. The error indicator is applied through the meshless RBF-FD method as found in the Medusa library [slak2021]. In general, the idea is to apply higher order explicit differential operators approximations to the implicitly obtained solution and thus indicate the areas with high error of the numerical solution. In the continuation of this work, the proposed error indicator will be named IMEX (implicit-explicit) error indicator.
Ii IMEX error indicator
Let there be a partial differential equation of type:
where is an arbitrary partial differential operator applied to , and equaling the constant . Such a problem is first solved implicitly, using a lower-order approximation of , , obtaining the solution in the process. The is then used to reconstruct explicitly with the help of higher-order approximation of , , giving . Finally, is then tested against the analytical to indicate the error. These steps can be summarized as follows:
compute approximations and ;
solve implicitly, obtain ;
compare and to indicate high error areas.
Iii RBF-FD approximation of differential operators
Since the introduction of meshless methods in the 1970s, many variants have been proposed. The first mention of RBF-FD dates from 2000 with the introduction from Tolstykh [tolstykh2000using]. Since then, the method has been thoroughly studied and applied to many real-world problems with recent applications to fluid flow [rot] and plasticity [filip] problems.
In the framework of RBF-FD, a linear differential operator in the node is approximated over a set of neighbouring (often called stencil) nodes
for an arbitrary function and weights yet to be determined. The weights are obtained by constructing a localised RBF approximation with a given set of radial basis functions (RBFs) centred at the stencil nodes of a central node
The localized intepolation (2) can be written in a linear system
However, as previously observed by Bayona et al. [bayona2017role], RBFs alone do not guarantee convergent behaviour or solvability of the system. To mitigate these problems, the approximation basis is extended by a set of monomials with up to and including degree in a -dimensional domain.
With the additional constraints, the RBF-FD approximation can be written compactly as
where is a matrix of monomials evaluated at stencil points,
is the vector of values composed by applying the operator under considerationto the polynomials at , i.e. and are Lagrangian multipliers (which we discard after the solution had been obtained).
The IMEX error indicator’s performance is demonstrated on a synthetic example, which is commonly used when testing adaptive algorithms in mesh-based methods [mitchell2013].
The example is the Poisson equation, which is solved in a 2D circular domain with its center at (0, 0), and radius 1:
The Neumann, and Dirichlet boundary conditions are defined through , and , respectively:
From these one can derive the analytical solution of the Laplacian at point :
The source is positioned at while controls the source strength. is the boundary normal at on . is positioned at (0.5, 0.5), and is set to 1000.
The example was solved on a laptop with Intel Core i7-8750H CPU, and 16 GB RAM. Results were computed, and written into a file in about 2 s111The source code for the example can be found at: https://gitlab.com/e62Lab/2022_CP_splitech_IMEX_error_indicator_poisson_eg.
V Results and Discussion
The computational domain is discretized and filled with scattered nodes using Medusa’s built-in algorithms [slak2019, slak2021]. This procedure results in a domain discretized with 24882 points. An example solution is shown in Fig. 1.
Support sizes for , and are set to (following the recommendations by Bayona et al. [bayona2017role]), being the monomial degree, and the number of dimensions of the domain.
The system in Eq. (6) is first solved implicitly, with lower order approximation of differential operators , and , which were obtained with 2nd degree monomials.
The solution for the scalar field is obtained with Eigen’s
BiCGSTAB solver [eigenweb].
To compute the RHS explicitly, a higher order approximation of the operator , obtained with 4th degree monomials, is applied to .
The results are then compared to produce the IMEX error indicator :
For validation purposes, the error of , , is also computed by comparing the implicit to the analytical solution. The latter is obtained with Eq. (8), and is:
. As the solution was obtained on scattered nodes, the source for the aforementioned line is obtained by Shepard interpolation (Python,
photutils.utils[larry_bradley_2020_4044744]), sampling each plot line point from 9 nearest neighbors.
Additionally, the same case is solved with 6th degree monomials used to produce for IMEX, with results plotted in Fig. 5.
Comparing Figs. 2, and 3 it is noticeable that the solution’s error is the biggest around the source at point . The IMEX error indicator also predicts the biggest error to be around the same point, as can be seen in Fig. 3. This is further supported by the graph in Fig. 4. Although the IMEX error indicator does not follow the actual error, it successfully identifies the area of the biggest error. Increasing the monomial degree to compute does not noticeably impact IMEX’s performance, as can be seen by comparing Fig. 4, and 5. However, increasing the monomial degree results in a significant compute performance hit in this particular case (total computation time increased to 4 s, compared to previous 2 s).
A synthetic example of the Poisson equation was solved and the IMEX error indicator was tested on it. The error indicator correctly indicated the area of increased error, which also coincided with the source in the Poisson equation. Results were produced by increasing the monomial degree of the explicit approximations by 2 compared to the implicit counterparts. Further increasing the monomial degree did not prove beneficial in this specific example.
We show that the proposed error indicator successfully identifies the areas with high error of the numerical solution. In the continuation, these findings could be used to adaptively refine the critical areas and improve the precision of the numerical solution.