I am aiming at developing a simulation software mostly aimed at reactive fluid flows applications. Here are a simplified list of points to drive the initial choices (code structure, libraries...):

  • I aim at implementing high-order finite-volume schemes, e.g. WENO5
  • I aim at using high-order time integrators, fully-implicit or IMEX or splitting
  • Ideally, the mesh adaptation should be driven by objective error estimates, e.g. multiresolution, so that spatial accuracy is properly ensured (instead of using other adaptation heuristic based, e.g. on gradients).
  • Native GPU capability would be an interesting bonus.

The recent library AMRex is an attractive choice, since it offers a comprehensive framework (mesh handling, time integrators, linear algebra) with GPU capabilities built-in as well. However, this library is focused on patch-based structured AMR. While this simplifies a great deal the mesh part of the implementation, this might be inefficient for cases where highly anisotropic refinement may be required (e.g. boundary layers).

My question is twofold:

  • What other considerations would you recommend adding to the list ? What libraries or framework would you suggest ?
  • In that particular boundary layer issue with complex geometries, what tricks do exist for structured adaptive meshes to obtain an anisotropic mesh in the physical domain near the desired boundary conditions ?

I am happy to detail my question more as the discussion unfolds. Thanks in advance for your help !

  • $\begingroup$ Reactive flow is a complex problem. Why would you want to write something yourself, rather than using what others have already spent dozens of person-years developing? Say, the PeleC class of codes from the US National Labs? $\endgroup$ Sep 19, 2023 at 14:10
  • $\begingroup$ There is added value in building your own solver, e.g. detailed control over the numerical schemes, the solution procedure, the implementation of application-specific models, etc... which motivate the creation of a specific platform in my case. $\endgroup$
    – Laurent90
    Sep 19, 2023 at 18:03
  • $\begingroup$ There is added value in building a demonstrator, but not in building another fully-fledged simulator. I will repeat that the Exascale Computing Project and its predecessors has probably spent several dozen person-years on building its reactive flow codes, which equates to several ten million dollars. No matter how much time you put into it, and no matter how talented you are, you cannot produce anything that is comparable. As others have pointed out below, you will get stuck on small pieces of the problem for a year each, and you will have nothing of interest to show at the end. $\endgroup$ Sep 20, 2023 at 2:51
  • $\begingroup$ I guess it depends on what "fully-fledged" means. I understand (but don't agree with) your point. Of course, building a complete simulation framework is a large undertaking, it is however perfectly feasible to develop specialised solvers for particular needs. Otherwise, why would there be so many solvers discussed in the literature, which seem to be more than just simple demonstrators ? Anyway, the motivation of my question is not the question itself ;) $\endgroup$
    – Laurent90
    Sep 20, 2023 at 9:39

1 Answer 1


Let me preface by saying you may have set on a path more adventurous than you had hoped... Mesh generation itself is often a bottleneck, let alone rigorous 4D adaptation of complex geometries powered by GPU. I have written in bold a few keywords which may prove helpful to search on scholar.

From my limited knowledge, here are a few (free) mesh adaptation libraries/pieces of software:

  • Inria's MMG (2D/3D simplices).
  • NASA's Refine (2D/3D simplices).
  • Avro (2D/3D/4D simplices & possibly quads/hexaedra).

These three are of the metric-based variety. A metric field (solution field of spd matrices) prescribes sizes and orientations for elements locally. Each of these remeshers then constructs a uniform mesh, in a space distorted by these metric fields.

These metric fields are derived by error estimators. The most basic rigorous error estimate controls the $L^p$ norm of interpolation error, these are also called multi-scale estimates (because of $L^p$ as opposed to $L^\infty$ capturing only largest scale variations), and part of sometimes called "Hessian-based mesh adaptation". This is because the driving error term (on linear meshes) of Lagrange $P^1$ interpolation error depends only on the Hessian of the considered solution field (which, in practice, is reconstructed numerically). See for instance A decade of progress on anisotropic mesh adaptation for computational fluid dynamics and references therein.

More advanced error estimators called goal-oriented or output-based control a functional of interest (say drag or lift) see e.g. here and here, or rely on local solves to provide generic a posteriori error estimates (which could prove more robust than a priori interpolation error estimates in the case of hyperbolic equations). See for instance here for results using these techniques (a posteriori goal oriented adaptation for CFD).

I believe at least MMG and Refine are capable of computing an optimal metric field for interpolation error directly from solution fields. For vector fields, I believe the standard is to compute a metric per component and to apply metric intersection.

If I understand correctly, you want space-time adaptation. Avro may be best suited to this, as it handles 4D meshes. However, I believe it only handles faceted geometries, and I don't know whether it implements any space-time error estimation, or if that's up to the user. I'm not very optimistic as to the fact another free piece of software will provide you with this capability either. Perhaps you will find this relatively simple method for interpolation error control or this one for an output error control reference helpful in implementing your own if you fail to find a ready solution. In these cases, space-time adaptation is done by, morally, adapting in space at each time step. Some "non-linearity" may be present, in that iteration cost targets may take into account the time evolution of the problem.

If you can opt for a non-conforming numerical scheme on quads, you can instead consider AMR libraries such as p4est or SAMURAI. To represent the geometry, you could either (least geometrically accurate solution) cut out elements intersection/outside the domain or (perhaps better but more involved) consider immersed boundary methods. In the latter, you compute explicitely the intersections between mesh elements and the boundary and apply a special treatment to those elements (i.e. modify the scheme).

A note regarding AMR; if the idea seduces you but you do not wish to implement a scheme on non-conforming meshes, and think hexaedra (or even quadrilaterals) may prove a decent comprise, I urge you to reconsider. For all their apparent simplicity, making meshes out of these elements is even more of a nightmare than out of triangles or tetrahedra.

Regarding boundary-layers specifically, you will most likely see them appear naturally using any of the above error estimator types, as evidenced by the numerical results in the linked papers. Still, it is not uncommon to be able to specify boundary-layer regions, in which case the remesher will generate them using frontal methods resulting in nice quasi-structured boundary-layers.

  • 1
    $\begingroup$ Thank you for the detailed answer. When I was mentioning time adaptation, I am only interested in dynamic on-the-fly time adaptation, without generating a time mesh per se. Simply said, the space discretisation (implemented with e.g. AMRex) yields a system of ODEs which I then then integrate forward in time. I'll take a look at the references you mentioned. $\endgroup$
    – Laurent90
    Sep 19, 2023 at 10:38

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.