Translate from libMesh meshes to MFEM meshes

[cmacmackin]

Closes #32520

[cmacmackin]

This is largely done. There are a few things that maybe could be tidied up in the code.

I'm also not totally satisfied with how the integration tests have been written. Ideally what I would like to do is run the same simulation twice: once with the mesh read directly in by MFEM and once with it first read by MOOSE/LibMesh and then converted to MFEM. The output Exodus output of those two simulations should be the same. Unfortunately, there didn't appear to be an easy way to do that. I've compromised by running one of these simulations as a parent app and the other in a sub-app and then finding the L2 error between the results. The problem with this is that, while it confirms the solutions are the same, it doesn't check that the meshes are. An alternative would be to save the result of loading the mesh directly into MFEM as a "gold" result, against which to compare the results when the meshes are converted. This is slightly more fragile however.

[cmacmackin]

I don't understand the nature of the current error when compiling in CI:

Compiling C++ (in opt mode) /tmp/build/moose/framework/build/unity_src/mfem_Unity.C...
In file included from /tmp/build/moose/framework/build/unity_src/mfem_Unity.C:134:
In file included from /tmp/build/moose/framework/src/mfem/utils/CubitElementInfo.C:1:
/tmp/build/moose/framework/build/header_symlinks/CubitElementInfo.h:2:10: fatal error: 'mfem/mesh/element.hpp' file not found
    2 | #include <mfem/mesh/element.hpp>
      |          ^~~~~~~~~~~~~~~~~~~~~~~
1 error generated.
make: *** [/tmp/build/moose/framework/build.mk:166: /tmp/build/moose/framework/build/unity_src/mfem_Unity.x86_64-pc-linux-gnu.opt.lo] Error 1
make: *** Waiting for unfinished jobs....

There seems to be a missing header file, but this is found just fine when compiling on my machine. I'm using the version of MFEM in framework/contrib.

Edit: I was forgetting to wrap some files in #ifdef MOOSE_MFEM_ENABLED. That fixed it.

[moosebuild]

Job Documentation, step Docs: sync website on 8ea808d wanted to post the following:

View the site here

This comment will be updated on new commits.

[cmacmackin]

There were a few different types of failure in the last CI

1. Running the integration tests with CUDA. I think this was just a problem in my configurations, which I have fixed.
2. The unit tests for converting from TRI7 -> TRI6 and TET14 -> TET10. The issue here is that doing so raises a warning, which the test suite is promoting to an error. I was able to write the tests in such a way to avoid this happening (by temporarily turning off this setting) on my machine, so I'm not sure why it is happening in CI.
3. Problems when running the integration tests in parallel. These errors are strange and depend on the particular combination of flags used with the test runner. It looks like they're related to MFEM and libMesh partitioning meshes differently, but I still don't fully understand the problem.

Additionally, the unit tests are complaining of a memory leak. When I tried running Valgrind on these locally it didn't report any problems. The leaks only happen when I work with MeshGeneratorMesh objects and seem to involved libMesh::Parameters::Value objects.

Edited 7 April to reflect which problems have now been fixed.

[lindsayad]

On your local machine are you using 64 bit or 32 bit dof indices? And what is the size of unsigned long?

[lindsayad]

I have pushed your branch to my fork after rebasing. It has a commit on top of it by codex which I haven't personally reviewed yet which at least removes the reported leaks from unit testing. I haven't attempted to verify whether this fixes the failing tests across all architectures yet either

[lindsayad]

Ok I reviewed the codex commit (lindsayad/moose@a8401461aae) and I think it's good. Main points

• Store an enum class ElementCase in the gtest's testing::Values instead of ElementTestData which contains InputParameters which ~LibmeshMeshInit will complain about because when in gtest's long-term Values storage, they outlive LibMeshInit. This admittedly requires more code but I think it's more idiomatic in terms of a libMesh/MOOSE test program
• Make sure to store element connectivity as a std::vector<dof_id_type> instead of std::vector<unsigned int> because that's the parameter type of the ElementGenerator
• Use distinct app/mesh etc. for the EXPECT_THROWS portion of the CheckWarning test. I think it's better not to rely on the data structures to be in a viable state post-throw and just proceed with new instances after the throw check

[lindsayad]

Ok I just pushed one more commit to my fork that was probably the root of the CI problem ...

template <typename S, typename... Args>
void
mooseWarningStream(S & oss, Args &&... args)
{
  if (Moose::_warnings_are_errors)
    mooseError(std::forward<Args>(args)...);

  std::ostringstream ss;
  mooseStreamAll(ss, args...);
  std::string msg = mooseMsgFmt(ss.str(), "*** Warning ***", COLOR_YELLOW);
  if (Moose::_throw_on_warning)
    throw std::runtime_error(msg);

If Moose::_warnings_are_errors, then it doesn't matter what you've set for Moose::_throw_on_warning. I found that I was seeing unit tests pass when I gtest filtered just the CheckWarning tests, but failed when running the full unit test suite. Somewhere we must set Moose::_warnings_are_errors

[cmacmackin]

Use distinct app/mesh etc. for the EXPECT_THROWS portion of the CheckWarning test. I think it's better not to rely on the data structures to be in a viable state post-throw and just proceed with new instances after the throw check

Wouldn't it just be simpler to split these into separate tests, so we can let the existing fixture handle this work?

[cmacmackin]

Ok I just pushed one more commit to my fork that was probably the root of the CI problem ...

template <typename S, typename... Args>
void
mooseWarningStream(S & oss, Args &&... args)
{
  if (Moose::_warnings_are_errors)
    mooseError(std::forward<Args>(args)...);

  std::ostringstream ss;
  mooseStreamAll(ss, args...);
  std::string msg = mooseMsgFmt(ss.str(), "*** Warning ***", COLOR_YELLOW);
  if (Moose::_throw_on_warning)
    throw std::runtime_error(msg);

If Moose::_warnings_are_errors, then it doesn't matter what you've set for Moose::_throw_on_warning. I found that I was seeing unit tests pass when I gtest filtered just the CheckWarning tests, but failed when running the full unit test suite. Somewhere we must set Moose::_warnings_are_errors

Huh. Very strange that this wasn't happening on my system, even when I ran the full test suite. I guess there must have been some part of the code I wasn't building that did this.

[loganharbour]

Use distinct app/mesh etc. for the EXPECT_THROWS portion of the CheckWarning test. I think it's better not to rely on the data structures to be in a viable state post-throw and just proceed with new instances after the throw check

Wouldn't it just be simpler to split these into separate tests, so we can let the existing fixture handle this work?

Separate tests or just enclose them in different scope within the same test. Either is fine

[cmacmackin]

But doesn't make parallel distributed work...

There was strange behaviour depending on what combination of parallel and distributed were used. Different combinations of tests failed in different ways. Most bizarrely, if you ran in serial but with the distributed flag (admittedly, a very odd thing to do) some of the tests failed.

[cmacmackin]

99072c2 makes parallel testing work

I'm not sure if this is actually the behaviour we want though. One of the reasons for this PR is so that we can implement efficient transfers between MFEM and libMesh meshes, therefore allowing a high level of interoperability between MFEM and libMesh based apps. Therefore, we were deliberately choosing to set up the MFEM mesh with the same partitioning as the libMesh one. Won't the changes in this PR mean that we can no longer guarantee the same nodes will be on the same ranks for both libMesh and MFEM meshes?

[lindsayad]

I think the premise of those multiapp copy transfer tests is fragile and worked even in the serial case by happy coincidence. They essentially rely on both libmesh and MFEM reading an exodus mesh and deciding on an internal element numbering representation that are the same. But libmesh or MFEM could have decided to take the exodus numbering and then reversed the numbering or done a totally random re-ordering and then there would be no chance of matching for copy transfers by just doing vector assignments

[lindsayad]

Given the test requirements, I don't understand why multiapps are involved at all

[lindsayad]

If we're going to do multiapp copy transfers the preparation process on parent and sub should be identical. Same initial Mesh type. Same communicator. Same distributed mesh flag if both originally libmesh

[lindsayad]

Keep in mind that I am no longer officially on this project so will have limited time to do any refactoring.

I'd be willing to finish up this PR if you want to pass it along

[lindsayad]

I'm not sure if this is actually the behaviour we want though. One of the reasons for this PR is so that we can implement efficient transfers between MFEM and libMesh meshes, therefore allowing a high level of interoperability between MFEM and libMesh based apps. Therefore, we were deliberately choosing to set up the MFEM mesh with the same partitioning as the libMesh one. Won't the changes in this PR mean that we can no longer guarantee the same nodes will be on the same ranks for both libMesh and MFEM meshes?

The only way in my opinion for you to generate a transfer as efficient as a copy transfer is if the meshes are indeed copies. And one app having type = FileMesh and another having type = MFEMMesh simply aren't copies unless the same exodus reading code translating to internal representation is the same, and they're not. So if the meshes aren't copies conceptually we just shouldn't be trying to use a copy transfer. You're going to have to use something less efficient in order to guarantee the transfer is robust

[alexanderianblair]

I think the premise of those multiapp copy transfer tests is fragile and worked even in the serial case by happy coincidence. They essentially rely on both libmesh and MFEM reading an exodus mesh and deciding on an internal element numbering representation that are the same. But libmesh or MFEM could have decided to take the exodus numbering and then reversed the numbering or done a totally random re-ordering and then there would be no chance of matching for copy transfers by just doing vector assignments

I agree with this (it's only because the MFEM exodus reader grew out of the Apollo version anyway that this happy coincidence exists...)

For any future MFEM<->libMesh CopyTransfer, we would need consistent partitioning between source and target; setting the MFEM partitioning to be the same as the libMesh partitioning of the mesh its built from seems cleaner to me than 99072c2 . But these CopyTransfers are beyond the scope of this PR; we fundamentally just want to test that the mesh as-built from MFEM directly or from an MFEM mesh built from a libMesh mesh are the same in the sense that they represent the same geometry/elements/order etc.

The two options for this PR that I can see are:  

1. We wait until Enable MultiApp transfers between MFEM and libMesh problems on CPU #32315 is in and swap out the MFEMCopyTransfers for the MFEMShapeEvaluationTransfers that would provide the 'less efficient, more robust' transfers between meshes of different partitioning

2. We compare ordered, serialised meshes built via both routes, which should render any need for any CopyTransfers moot.

From discussions with @cmacmackin , I'm in favour of option 2 (the idea for rewriting these tests alluded to above). This could be done using mfem::ParMesh::GetSerialMesh and mfem::Mesh::GetHilbertElementOrdering to produce the ordered serialised mesh, which could be written out as an MFEM mesh from both routes for comparison in a test.
This should result in stricter tests than option 1, and as a bonus, we would gain the ability to write out MFEM meshes, which would be useful for re-use of adaptively refined meshes.

[lindsayad]

For any future MFEM<->libMesh CopyTransfer, we would need consistent partitioning between source and target

It's even more fundamental of an issue than partitioning. There is no partitioning in serial. The internal element numbering can be completely different in theory between a libmesh mesh built from exodus vs. an mfem mesh built from exodus. You'd have to go through and do some geometric matching and then element renumbering to be able to trust that you have consistency between the two. I haven't actually checked to see whether that's done in the current code

[lindsayad]

This could be done using mfem::ParMesh::GetSerialMesh and mfem::Mesh::GetHilbertElementOrdering to produce the ordered serialised mesh, which could be written out as an MFEM mesh from both routes for comparison in a test.

I didn't see this when responding originally. This is definitely an elegant solution! I support option 2 leveraging that MFEM API. It seems like in that case we only need plumbing for telling the libmesh-based mesh to reorder since we already have the parameter reorder_mesh for MFEMMesh

[cmacmackin]

This could be done using mfem::ParMesh::GetSerialMesh and mfem::Mesh::GetHilbertElementOrdering to produce the ordered serialised mesh, which could be written out as an MFEM mesh from both routes for comparison in a test.

I didn't see this when responding originally. This is definitely an elegant solution! I support option 2 leveraging that MFEM API. It seems like in that case we only need plumbing for telling the libmesh-based mesh to reorder since we already have the parameter reorder_mesh for MFEMMesh

I'm most of the way done implementing this. Unfortunately, it's uncovered a couple of underlying problems with MFEM and/or libMesh.

1. Even after reordering, MFEM doesn't seem to reorder the record of the boundary elements.
2. libMesh and MFEM assume different choice of quadrature points when reading in higher-order meshes. MFEM assumes Gauss-Lobatto, while libMesh assumes uniformly spaced (Lagrange). This only becomes an issue at 3rd order, but it does affect EDGE4 elements.

I'm going to write my comparison script to filter out the parts of the mesh which we expect to differ due to these issues. It's not ideal but it's probably the best we can do at the moment.

[lindsayad]

1. Even after reordering, MFEM doesn't seem to reorder the record of the boundary elements.

This feels like something that should be an MFEM issue

2. libMesh and MFEM assume different choice of quadrature points when reading in higher-order meshes. MFEM assumes Gauss-Lobatto, while libMesh assumes uniformly spaced (Lagrange). This only becomes an issue at 3rd order, but it does affect EDGE4 elements.

EDGE4 is libMesh's only third order geometric element and I don't think anyone uses it, although I guess MFEM users are just the kind of people who would use it and consequently create EDGE4 elements in their libmesh mesh generators before converting to an MFEM mesh. This feels like a fixable issue. You could create a MOOSE issue and just error for now? This is an issue I'd be happy to tackle in a follow-on

[cmacmackin]

This feels like something that should be an MFEM issue

Agreed. For now we're just planning to ignore that section of the mesh file when making comparisons.

EDGE4 is libMesh's only third order geometric element and I don't think anyone uses it, although I guess MFEM users are just the kind of people who would use it and consequently create EDGE4 elements in their libmesh mesh generators before converting to an MFEM mesh. This feels like a fixable issue. You could create a MOOSE issue and just error for now? This is an issue I'd be happy to tackle in a follow-on

The bigger issue with this is we need to work out which set of quadrature points is actually the correct one to use. Does the Exodus documentation specify anywhere?

[lindsayad]

The nodes in the exodus mesh define the mesh geometry for both libMesh and MFEM. If you are using a Lagrange basis in libMesh, then that is fundamentally tied to the mesh nodes. For MFEM there is no link (I believe) between the interpolation nodes used to define solution variable finite element bases and the mesh geometry nodes. If you wanted to make the interpolation nodes for a variable evaluation match on an EDGE4 element with a third order Lagrange FE in libMesh, then I think you'd need to set a basis_type of Uniform for the MFEM finite element space

[lindsayad]

quadrature points are generally separate concepts from interpolation points/nodes both in MFEM and libMesh. You can choose to collocate them if you want if you want things like diagonal mass matrices

[moosebuild]

Job Precheck, step Clang format on 6dde692 wanted to post the following:

Your code requires style changes.

A patch was auto generated and copied here
You can directly apply the patch by running, in the top level of your repository:

curl -s https://mooseframework.inl.gov/docs/PRs/32566/clang_format/style.patch | git apply -v

Alternatively, with your repository up to date and in the top level of your repository:

git clang-format 6c5870f0e93f4b019b3d4a287f486d192217fefe

[cmacmackin]

I've rewritten the tests to do a direct comparison of meshes. Most of them work, but not those for tets and (when running in parallel) pyramids.

[loganharbour]

Pedantic request: the test names are long. the test already has libmesh_mfem_conversation in the name from the folder. We don’t need libmesh_mfem_conversion/LibMeshToMFEMMeshQuad for example. libmesh_mfem_conversion/quad should be plenty. And if it’s not, the folder name can be corrected to give more context. We should limit all tests in a folder anyway to a single feature.

Second thing: I’m not a fan of check.sh. Git diff won’t work with installed applications.

[alexanderianblair]

I've rewritten the tests to do a direct comparison of meshes. Most of them work, but not those for tets and (when running in parallel) pyramids.

I suspect this is related to the tet mesh orientation issues discussed in mfem/mfem#3847. I've found that

Suggested change

	mfem::Mesh serial_mesh = _pmesh.GetSerialMesh(save_rank);
	mfem::Mesh serial_mesh = _pmesh.GetSerialMesh(save_rank);
	serial_mesh.ReorientTetMesh();

appears to fix the failing tet mesh test failures. Unfortunately, this method has been deprecated for a while (since mfem/mfem#1046) - I haven't yet found a suitable alternative that gives the ordered serial mesh in this case as we require.

EDIT: Looks like this does get things working in serial, but the pyramid and tet tests in parallel are still failing

[lindsayad]

I've opened mfem/mfem#5302 for hopefully consistent boundary element ordering

[cmacmackin]

quadrature points are generally separate concepts from interpolation points/nodes both in MFEM and libMesh. You can choose to collocate them if you want if you want things like diagonal mass matrices

Sorry, I misspoke. I meant the basis points.

[cmacmackin]

Second thing: I’m not a fan of check.sh. Git diff won’t work with installed applications.

What do you mean by this? Would this be run somewhere without git? What would you recommend instead? Just plain-old diff?

[cmacmackin]

The nodes in the exodus mesh define the mesh geometry for both libMesh and MFEM. If you are using a Lagrange basis in libMesh, then that is fundamentally tied to the mesh nodes. For MFEM there is no link (I believe) between the interpolation nodes used to define solution variable finite element bases and the mesh geometry nodes. If you wanted to make the interpolation nodes for a variable evaluation match on an EDGE4 element with a third order Lagrange FE in libMesh, then I think you'd need to set a basis_type of Uniform for the MFEM finite element space

That is what I've done when converting from libmesh to MFEM. However, if MFEM reads the same file it would not use Uniform for the basis_type of the mesh geometry nodes. It appears that libmesh and MFEM are interpreting the same input data differently. Which one is correct?

[moosebuild]

Job Test, step Results summary on 8ea808d wanted to post the following:

Framework test summary

Compared against fa12ff0 in job civet.inl.gov/job/3740302.

Added tests

Test	Time (s)	Memory (MB)
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshPyramid/run	0.71	67.00
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshTet/run	0.67	68.72
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshWedge/run	0.65	149.16
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshHex2/run	0.62	101.27
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshTet2/run	0.62	113.11
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshQuad/run	0.60	81.34
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshHex/run	0.60	85.60
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshWedge2/run	0.58	127.91
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshWedge2/verify	0.07	0.00
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshTet2/verify	0.07	0.00
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshWedge/verify	0.05	0.00
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshHex2/verify	0.04	0.00
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshQuad/verify	0.03	0.00
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshHex/verify	0.02	0.00
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshPyramid/verify	0.01	0.00
mfem/mesh/libmesh_mfem_conversion.LibMeshToMFEMMeshTet/verify	0.01	0.00

Modules test summary

Compared against fa12ff0 in job civet.inl.gov/job/3740302.

No added tests

Run time changes

Test	Base (s)	Head (s)	+/-	Base (MB)	Head (MB)
stochastic_tools/test:web_server_control.stochastic_control/batch_reset_multi	2.05	4.01	+96.00%	241.32	237.74

[lindsayad]

That is what I've done when converting from libmesh to MFEM. However, if MFEM reads the same file it would not use Uniform for the basis_type of the mesh geometry nodes. It appears that libmesh and MFEM are interpreting the same input data differently. Which one is correct?

It depends on whether the nodes, after transforming to reference space, are uniformly spaced or not. So in reality it's liked tied to the file writer. If libMesh wrote the file, then the nodes will certainly be uniformly spaced. If MFEM wrote the file, I'd expect them to correspond to the basis type of the FESpace used to define the MFEM mesh nodes.