Using libCEED in MOOSE

[lindsayad]

We have spent the past few months evaluating different means to enable GPU acceleration for MOOSE. @NamjaeChoi has done some really great work writing MOOSE kernels using static polymorphism where before we used run-time virtuals for computing quadrature-point residuals and Jacobians. These "GPU" kernels appear very similar to our traditional "CPU" kernels which would be ideal as our application developers will need to understand what changes they need to make to enable their applications to run on GPU. As a note, he has been writing the GPU kernels with the help of kokkos so that we can support various backends.

Currently he is evaluating our finite element shape functions on CPU using our traditional finite element backend libMesh. However, if we can avoid having to rewrite/abstract our shape function evaluation routines in libMesh, that would be ideal. So we are wondering whether it could make sense for us to use libCEED for basis function evaluation and then pull its CeedVector data into our kokkos kernels? I don't know if anyone has considered using libCEED in such a piecemeal way. We are constantly developing our migration plan. I'm happy to share more thoughts/ideas/info to expand the discussion or provide clarification (on what feels like it's been a ramble/stream-of-consciousness)

@MengnanLi91 FYI

[jrwrigh]

Currently he is evaluating our finite element shape functions on CPU using our traditional finite element backend libMesh.

Assuming you're referring to CeedBasisCreateTensorH1Lagrange, the basis transformation matrices are being created on CPU as well. I'm not sure why you'd want to calculate those matrices on GPU? They're tiny compared to the problem size.

Or do you mean trying to use CeedBasisApply on device?

[NamjaeChoi]

Computing shape functions for quadrature points on reference coordinate does not need to be on GPU, but the problem is their transformation to the physical coordinate of each element. Also, shape functions are not necessarily computed for quadrature points, but they can be computed for arbitrary points in MOOSE, so we were wondering if we could take advantage from libCEED capabilities.

Probably CeedBasisApply is closer to what we are looking for, but we would just like to get the physical shape functions, preferably with a device API that returns the desired values for a single element and is callable inside parallel regions, so that we can choose between storing them in a huge vector or get them on-the-fly during computation.

[roystgnr]

At least in libMesh we lose "same shape functions" as soon as we hit p=3, with "flipping"/"rotating" of basis functions on edges/faces between elements of arbitrary orientation. That's just a permutation (possibly with sign changes), though.

I don't actually know what the original design decision was based on; I know John and Ben had collaborated a lot with Wolfgang and it could have been as simple as copying a deal.II decision. When I took my own pass at optimization, it was for Cahn-Hilliard calculations on Hermite (and Bogner-Fox-Schmidt and whatever you call the 3D version) elements (with no permutations or sign changes) on a structured grid (so other than AMR shape functions and derivatives weren't changing between elements) and caching seemed like the obvious solution for my case.

I never got to the general case, but that might have been because I was too dismissive after looking into an optimization strategy or two that worked for bilinear formulations but not for nonlinear ones ... or because the codes I was trying to optimize next were aerothermochemistry, so nonlinear that the initial profiling showed something like half our time spent inside exp(), with the algebraic solve second and shape function derivatives not even on the map.

[lindsayad]

Is there an example of libCEED with DG?

[jeremylt]

If one exists, it would probably be in the MFEM repo?

[lindsayad]

where is that?

[jeremylt]

Woops, sticky keyboard. MFEM repo: https://github.com/mfem/mfem
Or maybe the Palace repo: https://github.com/awslabs/palace

I am not sure if they have such an example, but that's where it would possibly be

[lindsayad]

A couple reasons i feel for keeping native MOOSE q-functions (what we call kernels in MOOSE) as opposed to going whole-hog into libCEED:

1. we can design them in a way that feels familiar to our application developers (mentioned in the OP)
2. allows us to support full matrix assembly since many of our applications require a sparse matrix for their choice of preconditioning

I will demonstrate my complete ignorance and ask is it silly to do full matrix assembly on device just to copy the memory to host for the linear solve depending on the PC? In which case we might as well just focus on matrix-free for GPU and go 100% with libCEED? Or is the matrix-free focus of libCEED based more on the focus on high order finite element methods?

@NamjaeChoi talked about evaluating shapes based on inverse mapping from physical point location. We do that for DG methods when evaluating the shape functions for neighbor elements. That’s also central in our contact methods.

Does libCEED support trace space bases for e.g. HDG methods?

[jrwrigh]

What are the practical issues with going matrix free? Obviously a lot of preconditioners don’t work matrix free.

That's the biggest practical issue as far as I'm aware, though @jeremylt and @jedbrown may have more detailed thoughts.

Generally you'll have to do some amount of assembly for most preconditioners (p-multigrid not withstanding), though since it's just for preconditioning you can lag the assembly to amortize the cost. I've found that to work well even on CPUs in some cases; do matrix-free matvecs, but lagged full-matrix assembly for better preconditioning can net better performance than using full-matrix assembly for matvecs and preconditioning. That's problem dependent though.

What about for problems with advection character? ... What does HONEE do for preconditioning for Navier-Stokes? In practice I suppose that even for high Reynolds flows, whether doing e.g. RANS or implicit LES, you may have enough symmetric character that multigrid still works well…?

Currently we just do (lagged) variable point-block Jacobi preconditioning. libCEED allows for assembling just the point-block diagonal (so we never assembly the full matrix), which then becomes "variable" once strong BCs remove some DoFs from the global linear system. HONEE's current focus is on scale-resolving simulations, specifically DNS right now, and VPBJacobi works well enough for that. YMMV for RANS or LES. But we lag the preconditioner, often for multiple timesteps.

Obviously, if we run with an explicit timestepping scheme and use a lumped mass matrix, we don't need to worry about any preconditioning and then matrix-free flies. We're looking into that as well, though we're focused more on implicit timestepping at the moment.

As for p-multigrid, it's something we're curious to try out and see, but we don't think it'll have as large an impact as it does for solid mechanics. @jedbrown probably has more refined thoughts on that.

I had not heard of HONEE before this discussion.

I'd be surprised if you had. 😆 It's only been a separate project for a few months now, but was forked off of the fluids example in libCEED (hence why HONEE's commit history is older than a few months). There's a lot of architecture changes that I want to make to it before it's truly ready for others to jump in, but there are priorities with Aurora ESP for the code right now. Oh, and also me graduating. I guess that's a priority too. 😅.

Oh, and regarding HONEE and matrix-free performance, I forgot that I uploaded the SIAM CSE presentation with that data. It's here, see slide 10 and 11. Note that in order to get the breakdown in timestep cost, we had to turn on GPU timing which slows down the code a decent bit. Without that timing, we were seeing about 0.5 seconds per timestep for Q3 on that specific problem.

[jedbrown]

Indeed, p-MG is natural for solid mechanics, but preconditioners for fluids are rather regime-dependent. We've had pretty good luck for compressible exterior flows (including low-Mach) using VPBJacobi, but steady flows/RANS and some low-Mach/incompressible interior flows will require multigrid. There are a few approaches for creating reduced elliptic systems (e.g., [this](https://doi.org/10.1137/090775889 compares several variable formulations and associated reductions for atmospheric flows), either as part of semi-implicit time integration or block preconditioning of implicit or IMEX (like ARK or Rosenbrock-W). We're interested in making quantitative comparisons of such formulations, but have been focused on the scale-resolving regime first.

[lindsayad]

More on the subject of fluids, how do you handle saddle point problems like incompressible Navier-Stokes matrix-free? I've spent much of the last year working on Schur-complement based preconditioners like LSC and AL using -pc_type fieldsplit -pc_fieldsplit_type schur

[rezgarshakeri]

We have the fieldsplit setup for incompressible elasticity in Ratel given in this example. You can find other cases like block diagonal there (see yml files with pcfieldsplit). In upper triangle preconditioner case [A B; 0 D]^-1 we have two separate Qfunctions for displacement block A and pressure block D ( we don't compute schur complement in pressure block) where we use p-MG on the first one and jacobi on the second.

[jedbrown]

Yup, most algorithms for saddle problems can be implemented matrix-free so long as the submatrices that are returned can procure entries when needed (like in the p-MG/AMG on the momentum block and/or the pressure block).

[hugary1995]

I am starting a new sub-thread so as not to divert previous discussion. I am a user of the MOOSE solid mechanics module. The types of material we work with all involve some form of plasticity, and so I am most interested in what the p-mutligrid AMG preconditioner can bring to us for such problems. In terms of scalability, we've had the best luck with Hypre so far. The benchmark we've run is remarkably similar to the Schwarz Primitive example as described in the arXiv paper. The major difference is that our material is hyperelastic-plastic with linear elements. In our benchmark, the linear solve (with Hypre as the preconditioner) very quickly becomes the bottleneck. Therefore, I am interested in learning how the performance improvements (as illustrated in the arXiv paper) carry over to plasticity problems, or more generally how libCEED can help in our applications.

[jrwrigh]

I'm not sure about performance of p-multigrid with hyper-elasticity, but Ratel supports both. @jeremylt can comment more on that.

[rezgarshakeri]

Ratel covers linear plasticity, different hyperelastic materials (Neo-Hookean, Mooney-Rivlin, Ogden) in compressible and incompressible limits, phase-field modeling, ... and has a demo example for linear plasticity but we don't have performance study results in plasticity.

[lindsayad]

I presume if we are using a nonlinear solution method like Newton, that the matrix-free approach of libCEED is more efficient than the matrix-free approach where J or its action can be approximated via finite differencing residual evaluations? It certainly has the advantage of not being sensitive to a differencing parameter, which can be very difficult to select for multiphysics problems whose function may be more noisy for some physics than others

[jedbrown]

That's good intuition, though there are more factors

• Preconditioning: MFFD doesn't help you with a preconditioner, which you will almost always need in some sense (if only Jacobi).
• Often a tangent model can be cheaper. For example, the Jacobian for nonlinear diffusion can pull back the effective diffusivity tensor back to reference configuration so there is no geometry present when evaluating the Jacobian action. For hyperelasticity, the Jacobian action can be defined entirely in current configuration so you don't need to access or store any data about the initial configuration and deformation gradient. This isn't always true, and there are circumstances (especially if you use a naive implementation) when the Jacobian seems more expensive to apply than the nonlinear residual. In such cases, you might use a quasi-Newton method.

[jrwrigh]

Another factor is accuracy differences between an analytic Jacobian and MFFD. MFFD is limited to $\sqrt{\epsilon}$ (so 1e-8 for double precision) relative accuracy, which can be significant for stiff problems.

[lindsayad]

Definitely. And simultaneously a differencing parameter selection of 1e-8 with double precision is often too tight for some of our applications whose functions are noisy. I frequently have to set -mat_mffd_err 1e-5 or so for some of our apps to get more robust solves