I'm trying to learn about numerically solving PDE by myself.

I've been beginning with finite difference method(FDM) for some time because I heard that FDM is the fundament of numerous numerical methods for PDE. So far I've got some basic understanding for FDM and been able to write codes for some simple PDE lay in regular region with the materials I found in the library and Internet, but what's strange is, those materials I've got usually talks little about the treatment of irregular, curved, strange boundary, like this.

What's more, I've never seen a easy way to deal with the curved boundary. For example, the book Numerical Solution of Partial Differential Equations - An Introduction (Morton K., Mayers D), which contains the most detailed discussion (mainly in 3.4 from p71 and 6.4 from p199) I've seen until now, has turned to a extrapolation that is really cumbersome and frustrating for me.

So, as the title asked, as to the curved boundary, usually how do people deal with it when using FDM? In other words, what's the most popular treatment for it? Or it depends on the type of PDE?

Is there a (at least relatively) elegant and high-precision way to deal with the curved boundary? Or it's just an inevitable pain?

I even want to ask, do people actually use FDM for curved boundary nowadays? If not, what's the common method for it?

Any help would be appreciated.


3 Answers 3


Answering your last question first, do people actually use FDM for curved boundary nowadays I'd say the answer is no. In the commercial CFD world, 2nd order accurate finite volume schemes are the de-facto industry standard. One of the advantages of FV (and finite element/discontinuous galerkin approaches Jed mentioned) over FD is the much more natural handling of complex boundaries. FD does provide the foundation of a lot of numerical methods (FV included) and it's necessary to learn as a first step, but it's not advisable for large-scale complex problems.

As for dealing with complex boundaries in FD, I can think of two canonical ways, one of which is the interpolation/extrapolation method you mentioned. The other is to use body-fitted grid points in physical $(x,y)$ space with a conformal mapping to "computational" $\xi = \xi(x,y),\eta = \eta(x,y)$ space where $\Delta \xi = \Delta \eta = constant$. Then one can re-write terms like

$$\frac{\partial u}{\partial x}=\frac{\partial u}{\partial \xi} \frac{\partial \xi}{\partial x} + \frac{\partial u}{\partial \eta} \frac{\partial \eta}{\partial x}$$

where the $\frac{\partial (\xi,\eta)}{\partial (x,y)}$ terms are called metric terms and can be computed at the beginning of a problem (or for a simple domain you might have an exact conformal mapping available), and the $u$ derivatives can be computed on a logically simple computational domain. This process makes the implementation of boundary conditions straightforward, but it requires generating a sufficiently smooth, nominally orthogonal curvilinear mesh.

I'd say this body-fitted grid approach is the "most popular treatment" for dealing with curved boundaries in FD, with the caveat that FD methods themselves are not very "popular" anymore for complex applications. It's rare to see them still come up in CFD literature except for on very simple domains.

  • $\begingroup$ Your statement "I'd say the answer is no" is not correct. Visbal and Gaitonde work extensively with higher-order FD in the FDL3DI code. Also, NASA's OVERFLOW code is a FD code (as far as I know/can tell). $\endgroup$ Commented Dec 10, 2013 at 13:50
  • $\begingroup$ OVERFLOW was originally purely FD, but now it generally uses FV flux splitting (AUSM, HLLC, etc., in Ch 1 of your link.) It's also definitely "legacy" code. That FDL3DI link is from work in the 90s when high-order finite-element/DG based work was in its infancy and there weren't any demonstraby viable high-order-accurate finite volume schemes. I think you'd be hard-pressed to convince someone in 2013 to start development of a code based on that work's compact finite difference strategy. As elegant as it is, it's very restrictive for applications. $\endgroup$
    – Aurelius
    Commented Dec 10, 2013 at 14:12
  • $\begingroup$ I kind of disagree with the generality of your statement that it's not advisable to use FD for large-scale complex problems. Nowadays, people in HPC tend to recast their finite element schemes in a stencil-like fashion and use (semi-)structured grids to efficiently implement matrix-free solvers for extreme scale computing. Thus, as unfashionable as they are, people still actually want to use finite differences. Not to mention that there are applications where you can get away with structured meshes. For complex geometries standard FD is painful though and maybe that's what you wanted to state. $\endgroup$ Commented Dec 10, 2013 at 14:35
  • $\begingroup$ For simple curved geometries, high-order FD will win over high-order spectral-difference/volume, flux-reconstruction, or DG methods on an efficiency basis (accuracy/time). For complex ones, grid generation may be enough of a headache to make you try alternative approaches. One should not forget that the very considerable flexibility of the above-mentioned methods comes at a considerable cost, see this paper by Loehner. This is one reason why FDL3DI and OVERFLOW still see use. $\endgroup$ Commented Dec 10, 2013 at 15:08
  • $\begingroup$ @ChristianWaluga yes that's basically what I was trying to state. Obviously FD ideas find their way into other applications (e.g. gradients in FV being computed by finite differences), and in certain areas like DNS on simple geometries you see them used. But for general-purpose codes the trend over the last 2 decades has been pretty clear away from pure FD. $\endgroup$
    – Aurelius
    Commented Dec 10, 2013 at 15:28

Curved boundaries are covered in most CFD books, e.g., Chapter 11 of Wesseling or Chapter 8 of Ferziger and Peric.

While not a fundamental theoretical problem, the practical complexity of implementing boundary conditions for high-order methods on curved boundaries is a significant reason for interest in more geometrically-flexible methods such as the finite element method (including discontinuous Galerkin). Structured finite difference and finite volume grids are still used in some CFD simulations, but unstructured methods are gaining popularity and the local operations used by high-order unstructured methods are actually quite efficient, and thus may not suffer much loss in efficiency compared to similar FD methods. (Indeed, the geometric flexibility often makes them more efficient.)

  • $\begingroup$ Great answer Jed. There is a very step-by-step walkthrough of how to treat irregular BC's in a fluids problem found in my thesis p38-46. Frankly it is a major A*# pain to do this in FD formulations. The important insight to take is that curved BC's can be approximated by a large number of infinitesimal straight ones. $\endgroup$
    – meawoppl
    Commented Dec 10, 2013 at 17:40
  • $\begingroup$ @Jed Suppose that we are going to solve a CFD problem with irregular boundaries; Sticking to a $2^{nd}$ order (not high order) method both for an unstructured FVM and a structured FDM, can we still say that for the unstructured FVM implementation is easier and local operations are more efficient? $\endgroup$
    – Naghi
    Commented Jul 1, 2020 at 15:42

I have worked on high precision fdm for the past n years. and I have used electrostatics -2 dim laplace's equation as the example for explicitly developing the high precision algorithms. until about 4 years ago the problems were constructed with horizontal or vertical lines points of potential discontinuity. if you google my name and fdm high precision you should find the references. but this is not your question. your question is fdm and curved boundaries. about a year ago I presented an order 8 solution in hong kong(see A Finite Difference Method for Cylindrically Symmetric Electrostatics having Curvilinear Boundaries) which created order 8 algorithms for interior points close to the boundary and these would require of course points on the other side of the boundary. the points on the other side of the boundary were put there by simply extending the mesh to the other side. having done this the question was how do you find the values of these points when relaxing the mesh. it was accomplished by integrating from the boundary (known potential) to the point using the algorithms. it was reasonably successful and reasonably accurate ~<1e-11, BUT required 103 algorithms each individually crafted and it was somewhat brittle, unstable geometries could be found. to remedy the above a solution has been found(order 8 and below) using (one!) minimal algorithm and the solution exhibits a considerable robustness. it has been submitted but would be available as a preprint by emailing me. I believe this technique would be extensible to time independent pde's(linear required) other than laplace and to dimensions higher than 2. I have not considered the time dependent problem but the technique being a power series technique should be adaptable and applicable. david

  • 2
    $\begingroup$ If you could submit your paper to a preprint server (like arXiv, for instance), and then link to it here, that would improve your answer. Generally speaking, answers should not contain e-mail addresses. I also encourage you to make your answer more concise. $\endgroup$ Commented Dec 18, 2013 at 23:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.