Trouble with AMR/reinit with subdomain restricted variables in TransientExplicitSystem

[ajamar1]

• I am using libmesh 1.3.0, but I think I have been fighting this issue unknowingly for years.
• The trouble is with the processing of old and older solution vectors in a TransientExplicitSystem when calling EquationSystems::reinit() after a call to MeshRefinement::refine_and_coarsen_elements() for parallel problems only.
• Further, the issue is isolated to when I have subdomain restricted variables. For example, subdomain 1 has variable 0, and subdomain 2 has variable 1. The two subdomains interface such that common nodes should have DoFs associated with both of the two variables.
• After a call to AMR/reinit, part of the processor boundary changes to be aligned at the interface of the two subdomains.
• When all is said and done, the data in the old and older solution vectors is corrupted (value set to zero) at the nodes on that specific part of the processor boundary. The data in the current solution looks ok.

*This goes away if I run on 1 processor.
*This goes away if I put both variables in the entire domain
*I am seeing this on a TransientExplicitSystem, and I am not sure if it would happen on a TransientImplicitSystem.
*I am not convinced that its an issue associated with AMR. Rather, I think it just took AMR to move the processor boundary around such that is aligned at the variable interface. So I think its an issue embedded in all the reinit logic.

I have attached a stand-alone code, which manifests the bug, in the attached tarball. The tarball includes a README file, Makefiles, an example problem with tecplot layouts to view the important output files, and an archived results folder to show you what I was getting on my system.

Inevitably, I am sure you will be asking for more details. Roy knows how to reach out to me, and I am happy to discuss over the phone. I hope I am doing something wrong and that its not a bug.

libmesh_bug_report.tar.zip

[roystgnr]

Thank you very much; this is enough for me to reproduce the problem.

I'm afraid I'm struggling trying to figure it out, though!

I am using libmesh 1.3.0, but I think I have been fighting this issue unknowingly for years.

I can replicate it with the current master/HEAD.

I am not convinced that its an issue associated with AMR. Rather, I think it just took AMR to move the processor boundary around such that is aligned at the variable interface. So I think its an issue embedded in all the reinit logic.

This appears to be correct.

I hope I am doing something wrong and that its not a bug.

It took me a while to be sure (because I don't use TransientSystem and I'd forgotten about the weird unique_ptr-that-doesn't-actually-manage-memory hack we'd used to get it working with AMR while still backwards API compatible) but no, it doesn't look like you're doing anything wrong, it looks like a bug.

It really is only affecting the old and older solution, not the current solution, which baffles me, since those all go through the same projection code. Somehow the projection code has to be working on PARALLEL type vectors but not on GHOSTED type? I can't imagine how but that's the only lead I see at the moment.

[ajamar1]

Roy (and all), I appreciate the rapid support. Let me know if I can provide any additional information. I look forward to getting past this one. Adam On Oct 24, 2018, at 2:15 PM, roystgnr <notifications@github.com<mailto:notifications@github.com>> wrote: Thank you very much; this is enough for me to reproduce the problem. I'm afraid I'm struggling trying to figure it out, though! I am using libmesh 1.3.0, but I think I have been fighting this issue unknowingly for years. I can replicate it with the current master/HEAD. I am not convinced that its an issue associated with AMR. Rather, I think it just took AMR to move the processor boundary around such that is aligned at the variable interface. So I think its an issue embedded in all the reinit logic. This appears to be correct. I hope I am doing something wrong and that its not a bug. It took me a while to be sure (because I don't use TransientSystem and I'd forgotten about the weird unique_ptr-that-doesn't-actually-manage-memory hack we'd used to get it working with AMR while still backwards API compatible) but no, it doesn't look like you're doing anything wrong, it looks like a bug. It really is only affecting the old and older solution, not the current solution, which baffles me, since those all go through the same projection code. Somehow the projection code has to be working on PARALLEL type vectors but not on GHOSTED type? I can't imagine how but that's the only lead I see at the moment. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#1919 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AqW_pbln_00v45g1gPIyBifnr6dHhoGhks5uoLxIgaJpZM4X2wvu>.

[roystgnr]

Yeah, that might be it - after changing the library to use PARALLEL vectors for old/older solutions, I no longer see discrepancies between those solutions before-vs-after reinit.

And my "same projection code" claim was an overstatement - there's some localization in system_projection.C:262-315 which we might be getting right (or at worst inefficient) in the PARALLEL case but might be failing in your example for the GHOSTED case.

[roystgnr]

Let me know if I can provide any additional information.

Change your vector type in TransientSystem::add_old_vectors() from GHOSTED to PARALLEL and let me know what happens in your application code? You're probably actually using the ghosted degrees of freedom, in which case the answer will be "even worse corruption in opt mode, error in dbg mode", but if not then that change would be a good temporary workaround for you and a good data point for me. I'd hate to fix the problem your test case triggers only to discover there were actually multiple problems.

[ajamar1]

I don’t want to get ahead of myself, but are you declaring victory and suggesting I rewire something on my end or accept a patch from you? Or are you thinking a different path forward? Thanks, Adam On Oct 24, 2018, at 2:28 PM, roystgnr <notifications@github.com<mailto:notifications@github.com>> wrote: Yeah, that might be it - after changing the library to use PARALLEL vectors for old/older solutions, I no longer see discrepancies between those solutions before-vs-after reinit. And my "same projection code" claim was an overstatement - there's some localization in system_projection.C:262-315 which we might be getting right (or at worst inefficient) in the PARALLEL case but might be failing in your example for the GHOSTED case. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#1919 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AqW_pU7fOO2wJIosFwOksyPfVnw0AK8vks5uoL9CgaJpZM4X2wvu>.

[ajamar1]

Just got this. I will try it and let you know. Likely this evening. Adam On Oct 24, 2018, at 2:31 PM, roystgnr <notifications@github.com<mailto:notifications@github.com>> wrote: Let me know if I can provide any additional information. Change your vector type in TransientSystem::add_old_vectors() from GHOSTED to PARALLEL and let me know what happens in your application code? You're probably actually using the ghosted degrees of freedom, in which case the answer will be "even worse corruption in opt mode, error in dbg mode", but if not then that change would be a good temporary workaround for you and a good data point for me. I'd hate to fix the problem your test case triggers only to discover there were actually multiple problems. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#1919 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AqW_pUs96e50o8gMMgxoqFPd2sLsfTzUks5uoMAPgaJpZM4X2wvu>.

[ajamar1]

Roy, You were correct about things being worse if I am using the ghosted values. In my application code, I hit a seg fault when running in opt. Debug gives me a nice error No index 231 in ghosted vector. Vector contains [252,462) And empty ghost array. Adam On Oct 24, 2018, at 2:31 PM, roystgnr <notifications@github.com<mailto:notifications@github.com>> wrote: Let me know if I can provide any additional information. Change your vector type in TransientSystem::add_old_vectors() from GHOSTED to PARALLEL and let me know what happens in your application code? You're probably actually using the ghosted degrees of freedom, in which case the answer will be "even worse corruption in opt mode, error in dbg mode", but if not then that change would be a good temporary workaround for you and a good data point for me. I'd hate to fix the problem your test case triggers only to discover there were actually multiple problems. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#1919 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AqW_pUs96e50o8gMMgxoqFPd2sLsfTzUks5uoMAPgaJpZM4X2wvu>.

[roystgnr]

Just so you know I'm not vanished, and so I can warn you that this isn't going to be a quick fix, and so I can get my current thoughts in writing:

I can see multiple instances of the bug here (both in BuildProjectionList and GenericProjector) that would trigger in your case, I just haven't figured out how to fix them. It looks like the correct behavior of the system solution itself is a red herring - we're not being as efficient as we should there, but the inefficiency usually masks the bug, but I can't swear that "usually" is "always", so trying to get the same behavior from the older solution vectors isn't the right fix.

As I feared, the trouble is coming when one processor owns a node that supports a subdomain-restricted variable but doesn't own any elements that support that variable. Then when we do our local projections on the elements the processor does own, the restricted variable can get lost in the shuffle.

This isn't going to be fixable without a massive rearchitecture of GenericProjector. Fortunately:

1. Half the raison d'etre for GenericProjector was to be able to do such rearchitecture and get it automatically incorporated by all the use cases
2. I've had "rearchitect GenericProjector" on my todo list for a while, for performance reasons, so maybe we can kill two birds with one stone.

Instead of looping over elements first, then doing the vertices-then-edges-then-faces-then-interiors routine on each element, we can do that process in four loops; the extra loop iterator cost will be more than made up for by avoiding redundantly calculating the same node/edge/face data multiple times (e.g. 8/4/2 times for modal data on a hex mesh, and much more on tets!).

When we treat vertices as independent rather than as being attached to the element then we can look specifically at the variables supported on the vertex. This alone would fix the case you sent me, where there aren't any nodes newly created by refinement at the subdomain boundary.

But newly refined vertices are where we first run into trouble. It is entirely possible on a distributed mesh for a node to be owned by a processor that can't even see any elements in a subdomain which has subdomain-restricted variables on that node. Even on a replicated mesh the node owner may not have solution data on any of those elements. If that node is new and we can't just look up previous variable indices from an old_dof_object there, then we have to interpolate, and we have to look at an element in the subdomain to do so. So we're going to need a communication step here. Processors with nonlocal nodes on processor boundaries will want to check the node for subdomain-restricted variables, then if it finds those check the patch around the node to see if there are elements owned by the node owner and including all relevant subdomains, then if it doesn't find those do the interpolation itself and push the result to the node owner. (this will sometimes still be a false positive but it'll still be a lot less redundant than our current algorithm and hopefully asynchronous communication will avoid too much of that cost too)

We'll then have the same problem for edges, but the same new code ought to fix it. With faces we can "cheat" and skip communication because we always keep face neighbors ghosted (and we always will in cases where we're doing adaptivity and/or repartitioning).

[ajamar1]

Roy, Thanks for keeping after this and for the status report. It will take me a few reads to fully digest this. As we move forward, please let me know how I can help. I do not consider myself capable of altering libMesh code, but I can certainly help you test any potential fixes. Having this fixed will be very enabling for us. Thanks again, and have a good weekend. Adam On Oct 26, 2018, at 5:07 PM, roystgnr <notifications@github.com<mailto:notifications@github.com>> wrote: Just so you know I'm not vanished, and so I can warn you that this isn't going to be a quick fix, and so I can get my current thoughts in writing: I can see multiple instances of the bug here (both in BuildProjectionList and GenericProjector) that would trigger in your case, I just haven't figured out how to fix them. It looks like the correct behavior of the system solution itself is a red herring - we're not being as efficient as we should there, but the inefficiency usually masks the bug, but I can't swear that "usually" is "always", so trying to get the same behavior from the older solution vectors isn't the right fix. As I feared, the trouble is coming when one processor owns a node that supports a subdomain-restricted variable but doesn't own any elements that support that variable. Then when we do our local projections on the elements the processor does own, the restricted variable can get lost in the shuffle. This isn't going to be fixable without a massive rearchitecture of GenericProjector. Fortunately: 1. Half the raison d'etre for GenericProjector was to be able to do such rearchitecture and get it automatically incorporated by all the use cases 2. I've had "rearchitect GenericProjector" on my todo list for a while, for performance reasons, so maybe we can kill two birds with one stone. Instead of looping over elements first, then doing the vertices-then-edges-then-faces-then-interiors routine on each element, we can do that process in four loops; the extra loop iterator cost will be more than made up for by avoiding redundantly calculating the same node/edge/face data multiple times (e.g. 8/4/2 times for modal data on a hex mesh, and much more on tets!). When we treat vertices as independent rather than as being attached to the element then we can look specifically at the variables supported on the vertex. This alone would fix the case you sent me, where there aren't any nodes newly created by refinement at the subdomain boundary. But newly refined vertices are where we first run into trouble. It is entirely possible on a distributed mesh for a node to be owned by a processor that can't even see any elements in a subdomain which has subdomain-restricted variables on that node. Even on a replicated mesh the node owner may not have solution data on any of those elements. If that node is new and we can't just look up previous variable indices from an old_dof_object there, then we have to interpolate, and we have to look at an element in the subdomain to do so. So we're going to need a communication step here. Processors with nonlocal nodes on processor boundaries will want to check the node for subdomain-restricted variables, then if it finds those check the patch around the node to see if there are elements owned by the node owner and including all relevant subdomains, then if it doesn't find those do the interpolation itself and push the result to the node owner. (this will sometimes still be a false positive but it'll still be a lot less redundant than our current algorithm and hopefully asynchronous communication will avoid too much of that cost too) We'll then have the same problem for edges, but the same new code ought to fix it. With faces we can "cheat" and skip communication because we always keep face neighbors ghosted (and we always will in cases where we're doing adaptivity and/or repartitioning). — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#1919 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AqW_pU7Dw9PzbtmABPozsKcozfpOmn6aks5uo4eJgaJpZM4X2wvu>.

[ajamar1]

Roy, It looks like the correct behavior of the system solution itself is a red herring - we're not being as efficient as we should there, but the inefficiency usually masks the bug, but I can't swear that "usually" is "always", so trying to get the same behavior from the older solution vectors isn't the right fix Are you suggesting that the current solution may have a bug too? I was thinking about hacks on my end in an attempt to save you work. While it would be quite intrusive, I could put a wrapper class around 3 ExplicitSystems (not transient) and just treat them as current, old, and older while in reality they are all just current. I use several TransientExplicitSystems, so I’d have to change them all, but I think all I do with them is setting and getting, so I think I’d be safe so long as the bug doesn’t impact the current solution vector. I recognize that ultimately, you’ll want to fix libMesh, but I also recognize that you may have other priorities that will keep you from fixing this. Is this a valid approach, and would you recommend me going down this path, OR do you want to keep hacking away at the real problem? Thanks, Adam On Oct 26, 2018, at 5:15 PM, Amar, Adam J. (JSC-EG311) <adam.j.amar@nasa.gov<mailto:adam.j.amar@nasa.gov>> wrote: Roy, Thanks for keeping after this and for the status report. It will take me a few reads to fully digest this. As we move forward, please let me know how I can help. I do not consider myself capable of altering libMesh code, but I can certainly help you test any potential fixes. Having this fixed will be very enabling for us. Thanks again, and have a good weekend. Adam On Oct 26, 2018, at 5:07 PM, roystgnr <notifications@github.com<mailto:notifications@github.com>> wrote: Just so you know I'm not vanished, and so I can warn you that this isn't going to be a quick fix, and so I can get my current thoughts in writing: I can see multiple instances of the bug here (both in BuildProjectionList and GenericProjector) that would trigger in your case, I just haven't figured out how to fix them. It looks like the correct behavior of the system solution itself is a red herring - we're not being as efficient as we should there, but the inefficiency usually masks the bug, but I can't swear that "usually" is "always", so trying to get the same behavior from the older solution vectors isn't the right fix. As I feared, the trouble is coming when one processor owns a node that supports a subdomain-restricted variable but doesn't own any elements that support that variable. Then when we do our local projections on the elements the processor does own, the restricted variable can get lost in the shuffle. This isn't going to be fixable without a massive rearchitecture of GenericProjector. Fortunately: 1. Half the raison d'etre for GenericProjector was to be able to do such rearchitecture and get it automatically incorporated by all the use cases 2. I've had "rearchitect GenericProjector" on my todo list for a while, for performance reasons, so maybe we can kill two birds with one stone. Instead of looping over elements first, then doing the vertices-then-edges-then-faces-then-interiors routine on each element, we can do that process in four loops; the extra loop iterator cost will be more than made up for by avoiding redundantly calculating the same node/edge/face data multiple times (e.g. 8/4/2 times for modal data on a hex mesh, and much more on tets!). When we treat vertices as independent rather than as being attached to the element then we can look specifically at the variables supported on the vertex. This alone would fix the case you sent me, where there aren't any nodes newly created by refinement at the subdomain boundary. But newly refined vertices are where we first run into trouble. It is entirely possible on a distributed mesh for a node to be owned by a processor that can't even see any elements in a subdomain which has subdomain-restricted variables on that node. Even on a replicated mesh the node owner may not have solution data on any of those elements. If that node is new and we can't just look up previous variable indices from an old_dof_object there, then we have to interpolate, and we have to look at an element in the subdomain to do so. So we're going to need a communication step here. Processors with nonlocal nodes on processor boundaries will want to check the node for subdomain-restricted variables, then if it finds those check the patch around the node to see if there are elements owned by the node owner and including all relevant subdomains, then if it doesn't find those do the interpolation itself and push the result to the node owner. (this will sometimes still be a false positive but it'll still be a lot less redundant than our current algorithm and hopefully asynchronous communication will avoid too much of that cost too) We'll then have the same problem for edges, but the same new code ought to fix it. With faces we can "cheat" and skip communication because we always keep face neighbors ghosted (and we always will in cases where we're doing adaptivity and/or repartitioning). — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#1919 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AqW_pU7Dw9PzbtmABPozsKcozfpOmn6aks5uo4eJgaJpZM4X2wvu>.

[roystgnr]

Okay, on closer investigation:

It does appear to me that the gross-inefficiency-bug in PARALLEL-type vector projections is managing to cancel out the missing-ghosted-values bug, so that awful hack you propose ought to work. I'm working on a better solution but you might reasonably not want to wait that long.