Boundary corruption when loading ParMesh from files

[M8kmyday]

Boundary corruption when loading ParMesh from files

[M8kmyday]

Yes, it could be in ParSubMesh.

I can put together an example starting with a mesh, going to a ParMesh, going to a ParSubMesh, then ParPrint. That will be far more involved than my prior example and will take some time to code, debug, and test. I'll post back when I have that.

[M8kmyday]

I don't have a way to test vs. a true 2D ParMesh. All of my 2D meshes are generated by ParSubMesh.

[M8kmyday]

I put together the attached test package that demonstrates the issue starting with a mesh produced by gmesh. If you need to compile, there is a make file. If not, you can just go to the test directory and run "./run.sh". It sequences through 30 ranks and produces a pass/fail message for each run. The results I get are in the file "results.txt".
meshdemo.tar.gz

[sebastiangrimberg]

I have a likely culprit but I haven’t tested it. In Mesh::FinalizeTopology, there is a bug where elements are generated for a parallel mesh when a subdomain has no boundary elements but the parallel mesh does. This is probably triggered when you’re constructing your ParSubMesh. This bug is actually fixed as part of #2669, an old PR. You can cherry pick the changes involving ReduceInt in mesh.cpp and it should resolve your problem, I hope.

[najlkin]

Hmm, I see a different story here. I tried the code on 16 slots and the test truly fails, showing 11 boundary elements, but the files look totally fine to me with 8 boundary elements in total (also glvis sees that 🙂 ). I played with the test and it seems that the extra boundary elements are generated during finalization of the loaded mesh when refine or fix_orientation are true. @M8kmyday , could you confirm that with refine and fix_orientation being false the test passes on your machine?

[sebastiangrimberg]

You need to use the master branch of MFEM, exposing these arguments was introduced in #3918.

[M8kmyday]

Got it compiled, and I see the test setup working properly with no failures. I re-compiled my simulators and ran the 3D regression suite, and I do not see the failures I was seeing before, and nothing new cropped up. It looks like a 100% fix to me.

[najlkin]

Good to hear! 😉 Closing this one in favor of #4128 .

[najlkin]

Thanks for the demonstration, will have a look at it more thoroughly. However, I am not sure if it is a bug, but a feature 😅 ParMesh is a distributed mesh, which means that every peer has a piece of it. The methods for printing the mesh then do not reconstruct a monolithic mesh again, but simply collect all pieces of the mesh to a single file. This way, you may see how the mesh was split, which is sometimes very useful for debugging etc., but other times it may cause troubles...

[najlkin]

Check PrintAsSerial if it is not what you are looking for by chance 😉

[M8kmyday]

Ok, I'll look forward to hearing how it goes.

I need ParPrint because I use ParSubMesh to extract a 2D port mesh and write it to a file. Then the 3D program spawns a separate 2D simulator that imports the 2D mesh, solves, and writes out the computed port DOFs. Finally, the 3D program imports the DOFs to keep working. I need ParPrint to work consistently or I end up with re-ordered DOFs between the 2D and 3D programs that corrupts the 3D solution.

[najlkin]

So the question is, do you want to follow the partitioning of the 3D mesh in 2D or not? If I understand it correctly, you do not want to respect that in the 2D solver if it is parallel as well and you probably want to redistribute the whole thing among the processors. You are only concerned about returning from this 2D system back to 3D, right? If this is the case, try to output that "serial" mesh and solve the 2D problem with that. We will need to figure out then how to apply again the partitioning...

[M8kmyday]

Yes, I want the 2D solver to follow the 3D partitioning. I have to maintain the distribution of the mesh between the 2D and 3D solvers. Otherwise, the DOFs get re-ordered, and I cannot import the 2D solution into the 3D solution. If I use the serial mesh output, then the 2D solver will partition the mesh differently than the way the 3D solver did it, and the DOF orderings between the 2D and 3D solvers will be different. When I use ParPrint, then the imported mesh in the 2D solver distributes the same as the 3D solver, the DOF orderings remain the same, and my solution exported from the 2D solver reads into the 3D solver with the correct DOF order.

[najlkin]

I think you are making a bit too many deductions here, that DOF ordering is one thing and that is something we need to sort out, but more important is what you need for the code itself. As I read it, it does not make any sense for you to use the 3D partitioning in 2D, because the submesh will be distributed just among a small fraction of all processors (if you do not do some trickery with the partitioning). So again, please try to export the serial mesh from 3D and import it to 2D. If this works, and you get some reasonable result there, we may explore the ways how to get it back to 3D.

[M8kmyday]

I already went down the serial output path, and it doesn't work for the reasons that I explained. I exchanged emails with Mark Stowell on the issue, and he recommended that I switch over to ParPrint on Oct. 11, 2023. This methods works great except for the occasional corruption that I documented.

This is a fairly mature project. I presented my 3D work on OpenParEM3D at https://mfem.org/seminar/ on Jan. 9 about a month ago, and I presented my 2D work on OpenParEM2D at https://mfem.org/workshop22/ in Oct. 2022 after 1 year of work. The 2D simulator has an automated regression suite with 27 cases with something like 40,000 individual tests. The 3D simulator currently has an automated regression suite with 22 cases with about 750 individual tests.

I'm not looking to re-architect and re-write 30,000+ lines of code from a 2.5 year effort and about 5000 man-hours of labor to work around the ParPrint bug. Can we just get ParPrint fixed?

[najlkin]

It is not a bug of ParPrint, it saves the mesh as a parallel mesh, i.e. topologically discontinuous, which implies extra boundaries, nodes, faces etc. Maybe I am missing something here, but using the 3D partitioning in 2D is highly inefficient as I explained, so I am trying to help you here to improve it. Changing the way the mesh is loaded is not a big deal I believe, basically just swap the constructor of ParMesh. With that DOF reordering, I am willing to help and we may consider making some new methods for that.

Otherwise, if you insist to use principally incompatible partitioning in 2D then we need to find a way to identify the interior boundaries and you will have to filter them out in whatever procedures you do with them, which is not exactly a smooth solution.

One more thing to be sure we are on the same page, do you see those boundary elements already on the side of the 3D code before exporting the mesh?

[M8kmyday]

I'm not concerned about efficiency in the 2D solver since 2D solves represent a tiny fraction of the overall 3D solve times. However, more efficiency is better than less efficiency. I just don't want to spend a lot of time on it.

Per your question on boundaries, the ParSubMesh functionality creates the boundaries for the 2D meshes as needed. So yes, they are seen in the 3D code before exporting the 2D mesh.

My 2D program can read serial or parallel meshes. The issue is how to get the solved 2D fields into the 3D solver to apply as boundary conditions at the ports. The 2D tangential field solution must be forced at the driven port in the 3D space.

From 3D mesh to a 2D port mesh file:
import 3D mesh -> ParMesh -> ParSubMesh for 2D ports -> custom processing [rotate, flip (if necessary), set space dimension, drop z] -> save 2D mesh file

2D solver (frequency domain using Petsc/Slepc complex numbers):
import 2D mesh -> solve complex vector fields on the port -> Petsc complex Vec of DOFs -> Petsc VecView to file

From 2D solution to 3D boundary condition:
custom 2D import [DOFs] -> ParGridFunction (built off ParSubMesh fes) -> transfer to 3D ParGridFunction (built off ParMesh fes) -> 3D solve -> 3D ParGridFunction -> etc.

Going from 3D to 2D is no problem in either serial or parallel. The issue is how to get the 2D solution back into 3D. My understanding is that MFEM is built around the DOFs being in proper order in Vectors, so I have to get the ordering exactly right across distributed vectors. Right now, I rely on ParPrint to keep the ordering in place across common partitioning in the 2D and 3D simulators, but that works only almost all the time. Hence the point of this ticket. If there is another way to get the 2D solution back into the 3D space, then I can try that.

[najlkin]

I meant if you see those "non-boundary" elements in 3D, but I guess you do. That almost is then probably the case when the partitioning assigns the whole boundary to a single processor, so you have a serial mesh in 2D virtually. Am I right?

[sebastiangrimberg]

If I understand correctly, I think you just want to turn off ParMesh::print_shared (using ParMesh::SetPrintShared) in order to avoid printing the subdomain boundary elements along with the true boundary elements when printing the parallel mesh to disk. This is on by default because I think it is used by visualization programs like GLVis.

That said, it sounds like your entire workflow would be drastically simplified by just using SubMesh/ParSubMesh and the corresponding transfer methods for GridFunction objects, and not dealing with any file IO.

[M8kmyday]

I added a line for SetPrintShared:

     parSubMesh->SetPrintShared(false);
     parSubMesh->ParPrint(parout);

I see the same behavior as before.

[sebastiangrimberg]

If you just use ParPrint from a true (2D) ParMesh, do you get the same issue? Is this related to ParSubMesh or generally to parallel meshes?

[najlkin]

@sebastiangrimberg , that option SetPrintShared is only for Print, not for ParPrint. @M8kmyday , it seems to me that you have generally problems with discontinuity of the mesh for the 2D solver. If I am right that the code worked when the whole boundary was assigned to a single processor, which might be very often, then the simplest solution would to use some manually constructed/modified partitioning, which would guarantee this 😮

[M8kmyday]

My setup works fine when the 2D mesh is distributed across multiple processors. Different examples and setups distribute in different ways depending on how the mesh was split up with ParMesh and ParSubMesh. It is not the case that I sometimes get lucky and have all of the 2D mesh on a single processor. I get the same answer no matter how the mesh is split across processors in the 2D solution except when ParPrint makes interior mesh elements into boundary elements.

[M8kmyday]

In any event, my submission at the start of this ticket is very simple and has nothing to do with simulation. It shows three cases with the exact same mesh produced by ParSubMesh printed with ParPrint for 19, 20, and 21 slots. The case with 20 slots has interior elements marked as boundary elements. That is obviously going to cause problems no matter how a simulator is configured.

[najlkin]

Hmm, that breaks my theory 😄 . It might be a bug after all, but I do not see where it would come from in ParPrint, it must something with ParSubMesh itself perhaps. It is difficult to track it from just the output mesh files. Can you provide a piece of code producing the/some flawed meshes?

[najlkin]

Do not spend too much time on that, just something simple 😉

[jandrej]

I don't immediately recognize that this is caused by SubMesh, although it's possible. SubMesh and it's parallel version create a valid Mesh/ParMesh, the only additional information it stores is the mapping between itself and the parent.

Boundary elements in SubMesh are only created if the parent has boundary elements defined on those internal boundaries. If you create a 2D SubMesh from a 3D parent, all faces that are marked with a boundary attribute are going to be boundaries in the 2D mesh. Additionally, if the 2D SubMesh co-dim 1 element doesn't have an adjacent element it is considered a boundary element, this is necessary to guarantee a valid mesh.

How are those meshes created? Is the meshing utility inserting unwanted data?

[M8kmyday]

I don't think I understand your question. The 3D mesh is created in gmsh. After that, how the 2D mesh is created is documented in the source code of the example I provided. I'm not using any meshing utility outside of MFEM.

[najlkin]

@jandrej it seems it depends on the partitioning in MFEM and not on the initial mesh. I am looking into that...

[jandrej]

In Mesh::FinalizeTopology, there is a bug where elements are generated for a parallel mesh when a subdomain has no boundary elements but the parallel mesh does.

@sebastiangrimberg I think this is the intended behavior. You can't compute on the SubMesh otherwise.

[sebastiangrimberg]

You're right, it does look like ParSubMesh calls FinalizeTopology(false); and so the mesh partition boundary elements shouldn't be generated in any case. Perhaps though there is some issue with the ParSubMesh constructor creating extra 1D boundary elements in parallel when constructed from the boundary of a 3D mesh? Not sure about this, but that's where I suspect something is going on (should be easy enough for someone to trace this issue down at this point).

[najlkin]

Oh my, it is terrible 😩 . The thing is Mesh::GenerateBoundaryElements() is purely serial and does not check anyhow if the edge is shared or not. It is just luck if the loaded local mesh has at least one boundary element and then everything is ok and the original boundary elements are preserved even when the local mesh has many more shared edges. Otherwise when there is no original boundary element, the method generates a boundary element for every shared edge. This is totally inconsistent and not meaningful anyhow for distributed meshes 😔 .

[M8kmyday]

I was checking your request, but it does not look like it is needed now. Let me know if/when I can help.

[najlkin]

@M8kmyday , have you tried that refine and fix_orientation? If they are both false, that generation of boundary elements is skipped altogether. So please try it.

[M8kmyday]

I can't see how to set fix_orientation. There isn't a method, and it is not in a constructor. I can set it using ParMesh::Load, but for that I need to construct a ParMesh on a communicator, and I do not see how to do that. I can create a blank ParMesh, but I don't see a method to set its communicator. How can I set fix_orientation?

[najlkin]

It is in the constructor, try this ParMesh(PETSC_COMM_WORLD, ifs, false, 0, false);.

[M8kmyday]

It does not work for me:
meshtest.cpp:22:57: error: no matching function for call to ‘mfem::ParMesh::ParMesh(ompi_communicator_t*&, std::ifstream&, bool, int, bool)’

There are 5 arguments in your ParMesh call, but the longest ParMesh constructor only has 4. The constructor to load a mesh only has three arguments, and I can set refine there. That still leaves fix_orientation.

[M8kmyday]

I'm on the production release, and to switch to the Master branch will require a 1 hour re-compile, and then I won't be on the production release.

Since this is a near-trivial change to my test program, can one of you that are already on the Master branch make that change and try it?

[M8kmyday]

I'll go ahead and compile the Master branch. It can just grind away in the background while I work on other stuff.

[najlkin]

I have tried that of course (without PETSc though, just MPI) and it worked for me, so I wanted to check if it is the same problem you had or something different, because the description sounded somewhat different. However, it is not meant to be a solution for your main code, because you will probably experience problems with the mesh when orientation of the faces is not fixed, but you may try it. I have put together a solution, but I am testing it and wanted to check this with you first.

[M8kmyday]

make config is failing for me:

make[1]: Leaving directory '/home/briany/Desktop/mfem-4.6.2-rc0/config'
cd "." && ln -sf "/home/briany/Desktop/mfem-4.6/data" .
ln: ./data: cannot overwrite directory
make: *** [/home/briany/Desktop/mfem-4.6/makefile:647: build-config] Error 1

[najlkin]

You have some data directory there already...

[najlkin]

Ok, I just pushed the fix in #4128 . Try it (even with the old version of MFEM eventually) and let me know if it works. We would close this one in favor of the PR 😉 .

[M8kmyday]

I downloaded and compiled mfem-najlkin-pmesh-load-fix, but somehow, it is still pulling in parts of mfem-4.6. Can't seem to track down why. I'll try again tomorrow.