Lumped Mass Matrix Implementations

[jrwrigh]

Currently, all the lumped mass matrix implementations (in the fluids code at least) utilize a row-sum technique, which will not work for higher-order tets.

Possible solution to this would be to use "special lumping", as described in this thread: https://github.com/CEED/libCEED/pull/1147#discussion_r1111328658. Jed's comment sketches out a way to possibly implement that method at a global level rather than the element level presented in Hughes' FEM book.

Relevant existing threads:

• https://github.com/CEED/libCEED/pull/1147#discussion_r1111312341
• https://github.com/CEED/libCEED/pull/1149#discussion_r1146929182

[jedbrown]

A related question here is which of these projections need to be lumped. The diagonal is a fine preconditioner for a consistent L^2 projection, and typically has some quality improvements. I concede it may be too expensive for the SNES residual for the diffusive fluxes, and I think John mentioned at some point that lumping is arguably better there since it's more local.

I think the intent was to reference this comment for an implementation strategy that could go into PCJacobi.

[jrwrigh]

My general logic for whether lumped or consistent is whether the projection occurs every timestep, or is done less frequently. So the turbulent statistics is consistent (though that's also a very small problem to solve anyways since it's on a 2D slice), but the velocity gradient projection (needed for the data-driven SGS) is lumped since that's needed every non-linear iteration.

[jedbrown]

Let's try adding this conservative lumping to PCJacobi. Is that something you can do or should we recruit someone else? Should be just a few lines of code (MatMult, VecSum, VecScale).

[jrwrigh]

You mean as a built-in for PETSc, or building and setting the Pmat explicitly? I was originally thinking the latter, but I guess there's no reason that it can't be integrated into PETSc proper.

[jedbrown]

I'd make it -pc_jacobi_type conservative" or some similar name, alongside the current options diagonal,rowmax,rowsum. See include/petscpctypes.handsrc/ksp/pc/impls/jacobi/jacobi.c(andsrc/ksp/f90-mod/petscpc.h` to make it available in Fortran).

[knepley]

@jedbrown I cannot find a description of the lumping in that other PR. Is this something I can use for PyLith?

[jedbrown]

I don't trust it (because I haven't used it), but this is the page from Hughes' book (in collapsed comment just above previously linked one).

Here's the description of the method and comparison table in Hinton et al 1976.

[jrwrigh]

I have a functioning prototype that fails fluids regression tests by a "small", "acceptable" margin (at least for being a completely different method). Any specific tests you want to see?

[jedbrown]

Cool. You had some tests in a profile comparing consistent and lumped for hex elements. Can you do that with quadratic tets and include consistent, diagonal, and this lumping?

[jrwrigh]

You had some tests in a profile

Are you talking about the tests I did for the spanwise turbulent statistics collection? That should be easy to do for hexes, but the fluids code doesn't support tets right now.

[jedbrown]

Hmm, src/dm/dt/fe/tests/ex3.c incorrectly claims to do an "L2 projection", but is really just sampling the dual basis functions. @knepley what's your most convenient test of accuracy of a consistent mass projection (with a KSP solve that we can control as -ksp_type preony -pc_type jacobi -pc_jacobi_type conservative)?

[jrwrigh]

I tested an L2 projection on the blasius boundary layer for hexes. Q1 and Q2 have the same DoFs (so different element counts). I'm plotting the velocity profile normal to the wall, zoomed in to see the differences:

Q1, at wall (MeanVelocity_MeanVelocityX here is the conservative lumped mass matrix):

Q2, at wall:

Q2, mid boundary layer:

A few observations:

1. The rowsum and conservative lumped mass matrix appear to give more-or-less the exact same answers. Not sure if this is expected or not.
2. For Q2, the error at the wall goes away; we're able to recover the exact zero velocity at the wall for the lumped matrices.
3. For Q2, the lumped mass matrices away from the wall appear to almost be using the linear-form of the mass matrix; the resulting projection appears to be completely linear within the element bounds.

[jrwrigh]

I've confirmed that the conservative and rowsum entries are not identical to each other, but are still pretty close (order 1e-4 maximum relative difference).

[jedbrown]

Cool that they're that close. I don't a priori have an explanation for that. We could use the projection of diagnostic variables in Ratel as a test of interesting fields on tet meshes.

[jeremylt]

@jrwrigh is this line of effort over in HONEE now?

[jrwrigh]

It isn't, but should be. I'll make an issue over there and close this for now.

[jrwrigh]

Actually, I think the general idea was to implement it in PETSc, but it hasn't been a priority since we've only been working with tet elements.

[jedbrown]

I think you mean "only been working with hex elements". Note that with MatSetValuesCOO, it would be possible to do element-based lumping so long as it was computed before assembly (and we'd need to know element volumes).

[jrwrigh]

Ope, yes. We haven't been working with tet elements.

I had the thought too for the element-based lumping, but we can cross that bridge when we get there.