From: Tim A. <tra...@bu...> - 2016-11-21 16:43:57
|
Hello: Sorry for the longish email that follows, I tried to keep it as succinct as possible. tl;dr Is there any way to alter EquationSystems::reinit() to not call refine/coarsen_elements()? These calls prevent me from properly reinitializing my new QoI class. Or, is there another way to use the element refinement flags to identify elements that were just refined/coarsened? @pbauman pointed me to a GitHub PR <https://github.com/libMesh/libmesh/pull/1080> discussion where it seems that these refine/coarsen calls in ES::reinit() are not really necessary anymore. ----------------------------- I am currently working on implementing a new QoI in GRINS, but I am having an issue with EquationSystems::reinit(). My new QoI is based on the GRINS::RayfireMesh class, which must be reinitialized after every refine/coarsening. To fully integrate this, I am adding hooks into GRINS::MultiphysicsSystem::reinit() that will call reinit() on my QoI class. I have a test case that follows this flow: assemble QoI, refine certain elements, reinit the EquationSystems (which would also reinit my QoI), assemble QoI again. The reinit() <https://github.com/grinsfem/grins/blob/master/src/qoi/src/rayfire_mesh.C#L160> for RayfireMesh looks for an INACTIVE element with JUST_REFINED children in order to refine its own internal 1D mesh. The problem arises in that MeshRefinement::refine_elements() ends up getting called *twice*, once in my test case and again in ES::reinit() ( https://github.com/libMesh/libmesh/blob/master/src/systems/equation_systems.C#L265 ) When refine_elements() is called the second time, it replaces all JUST_REFINED flags with DO_NOTHING, breaking my RayfireMesh reinit function. If I do not call refine_elements() in my test, and instead just flag the elements and let the ES::reinit() do the refinement, then this <https://github.com/libMesh/libmesh/blob/master/src/systems/equation_systems.C#L270> IF statement is skipped, and prolong_vectors() below it trips an assert. Perhaps I could alter how I am identifying the just refined elements? The RayfireMesh map stores the IDs of main mesh elements along the 1D rayfire. On reinit(), it looks for inactive elements, and then tries to determine if they were refined or coarsened by looking at the refinement flags of its children/parent. Would there be any other situation where an element would go from ACTIVE to INACTIVE with DO_NOTHING children? I'm not sure how robust of a check that would be. I would appreciate any thoughts/insights/suggestions for how to address this. Thank you, -- *~Tim Adowski* |
From: Roy S. <roy...@ic...> - 2016-11-23 22:45:37
|
On Mon, 21 Nov 2016, Tim Adowski wrote: > Is there any way to alter EquationSystems::reinit() to not call > refine/coarsen_elements()? These calls prevent me from properly > reinitializing my new QoI class. Or, is there another way to use the > element refinement flags to identify elements that were just > refined/coarsened? > > @pbauman pointed me to a GitHub PR discussion where it seems that > these refine/coarsen calls in ES::reinit() are not really necessary > anymore. This requires some care. Let me think out loud. The current AMR in ES::reinit() is three steps: 1: for backwards compatibility, we assume that the user called for MeshRefinement themselves, and we run a prolong_vectors(). 2: coarsen_elements() (and a restrict_vectors(), if anything changed), to shrink the system as much as refinement flags allow. 3: refine_elements() (and a prolong_vectors(), if anything changed), to grow the system to its final requested size. In the long run, there's probably three classes of behavior we still want to support: A. Backwards compatibility with apps that call refine_and_coarsen_elements() themselves. B. Backwards compatibility with apps that flag elements themselves but expect libMesh to do the refinement and coarsening. We could do that two ways: B1. Two-step coarsen-then-refine, for apps that worry about maximum memory use. B2. One-step coarsen-with-refine, for apps that worry about total CPU use. And in all those cases, we want: α. The final refinement flags should properly indicate JUST_REFINED/JUST_COARSENED elements. You need us to support α. (α+A, if you're using the adaptivity code in GRINS, right?) We need to not break A or B. We can support either B1 or B2 in the short run without adding an API for both. When we're in coarsen_elements() or refine_elements(), we "clean up" JUST_REFINED and JUST_COARSENED flags. If we don't do that, and someone tries to use case B in the code, then they'll be in our B1 implementation, the projection code will be confused by the leftover flags, and we'll die. But if we want to get α+B1 to work at all, then we have to handle exactly this problem: when the refine step occurs, it *can't* clear the JUST_COARSENED flags. Perhaps we can give the project_vector code an option to ignore JUST_COARSENED or JUST_REFINED flags, depending on whether it's being called from restrict_vectors() or from prolong_vectors()? Then coarsen_elements will only clear old JUST_COARSENED flags (before it sets new up to date flags itself), refine_elements will only clear old JUST_REFINED flags (before setting new up to date flags itself), restrict_vectors() will ignore the JUST_REFINED case, prolong_vectors() will ignore the JUST_COARSENED case, and everyone will be happy. This isn't trivial, though, and it affects code paths which we need to support yet may have no test coverage of, and I'm behind on other things. If you're up to adding that option yourself, that would be fantastic, and I'll help with review and debugging. If you can wait a month and a half then that's okay, and I can get things working by myself. The best application-level workaround I can think of is: 1. Make sure you configured libMesh with unique_ids on. 2. Move MeshRefinement::_smooth_flags() to a public MeshRefinement::smooth_flags() function. 3. Call smoothing manually before refinement. 4. Save flags manually. I'd suggest a hash map from element unique_id to RefinementFlag. After the smoothing is done this should be easy: if an element is marked to coarsen then save JUST_COARSENED on its parent, if marked to refine then save JUST_REFINED on it (which isn't quite right, but the children don't exist yet so you don't know what their unique ids will be). Maybe there's an intermediate-difficulty option? It would be reasonable as a temporary workaround to add a loop over elements, in each of refine_elements() and coarsen_elements(), which tests for REFINE or COARSEN flags, and then *only* clears JUST_REFINED/JUST_COARSENED flags if there is new refinement or coarsening to be done. That would be a much easier library modification and it would fix the α+A use case, even though it wouldn't be a step in the right direction long-term. If you do this, though, you definitely need to add a unit test, so I don't accidentally break your code again when we ever do get around to long term refactoring. --- Roy |
From: Tim A. <tra...@bu...> - 2016-11-26 15:58:26
|
Thanks for the reply Roy. I figured this wouldn't be a trivial thing to address. I have an AIAA paper and a thesis to get done, so I think for now I will just try hacking my way around this to get my GRINS code working. If/when you do get to looking at this, I will have a few test cases that might help you. Thanks again! ~Tim On Wed, Nov 23, 2016 at 5:45 PM, Roy Stogner <roy...@ic...> wrote: > > On Mon, 21 Nov 2016, Tim Adowski wrote: > > Is there any way to alter EquationSystems::reinit() to not call >> refine/coarsen_elements()? These calls prevent me from properly >> reinitializing my new QoI class. Or, is there another way to use the >> element refinement flags to identify elements that were just >> refined/coarsened? >> >> @pbauman pointed me to a GitHub PR discussion where it seems that >> these refine/coarsen calls in ES::reinit() are not really necessary >> anymore. >> > > This requires some care. > > > Let me think out loud. > > > The current AMR in ES::reinit() is three steps: > > 1: for backwards compatibility, we assume that the user called for > MeshRefinement themselves, and we run a prolong_vectors(). > > 2: coarsen_elements() (and a restrict_vectors(), if anything > changed), to shrink the system as much as refinement flags allow. > > 3: refine_elements() (and a prolong_vectors(), if anything changed), > to grow the system to its final requested size. > > > In the long run, there's probably three classes of behavior we still > want to support: > > A. Backwards compatibility with apps that call > refine_and_coarsen_elements() themselves. > > B. Backwards compatibility with apps that flag elements themselves > but expect libMesh to do the refinement and coarsening. We could do > that two ways: > > B1. Two-step coarsen-then-refine, for apps that worry about > maximum memory use. > > B2. One-step coarsen-with-refine, for apps that worry about total > CPU use. > > > And in all those cases, we want: > > α. The final refinement flags should properly indicate > JUST_REFINED/JUST_COARSENED elements. > > > You need us to support α. (α+A, if you're using the adaptivity code in > GRINS, right?) We need to not break A or B. We can support either B1 > or B2 in the short run without adding an API for both. > > > When we're in coarsen_elements() or refine_elements(), we "clean up" > JUST_REFINED and JUST_COARSENED flags. If we don't do that, and > someone tries to use case B in the code, then they'll be in our B1 > implementation, the projection code will be confused by the leftover > flags, and we'll die. > > But if we want to get α+B1 to work at all, then we have to handle > exactly this problem: when the refine step occurs, it *can't* clear > the JUST_COARSENED flags. > > Perhaps we can give the project_vector code an option to ignore > JUST_COARSENED or JUST_REFINED flags, depending on whether it's being > called from restrict_vectors() or from prolong_vectors()? Then > coarsen_elements will only clear old JUST_COARSENED flags (before it > sets new up to date flags itself), refine_elements will only clear old > JUST_REFINED flags (before setting new up to date flags itself), > restrict_vectors() will ignore the JUST_REFINED case, > prolong_vectors() will ignore the JUST_COARSENED case, and everyone > will be happy. > > This isn't trivial, though, and it affects code paths which we need to > support yet may have no test coverage of, and I'm behind on other > things. If you're up to adding that option yourself, that would be > fantastic, and I'll help with review and debugging. If you can wait a > month and a half then that's okay, and I can get things working by > myself. > > > The best application-level workaround I can think of is: > > 1. Make sure you configured libMesh with unique_ids on. > > 2. Move MeshRefinement::_smooth_flags() to a public > MeshRefinement::smooth_flags() function. > > 3. Call smoothing manually before refinement. > > 4. Save flags manually. I'd suggest a hash map from element > unique_id to RefinementFlag. After the smoothing is done this should > be easy: if an element is marked to coarsen then save JUST_COARSENED > on its parent, if marked to refine then save JUST_REFINED on it (which > isn't quite right, but the children don't exist yet so you don't know > what their unique ids will be). > > > Maybe there's an intermediate-difficulty option? It would be > reasonable as a temporary workaround to add a loop over elements, in > each of refine_elements() and coarsen_elements(), which tests for > REFINE or COARSEN flags, and then *only* clears > JUST_REFINED/JUST_COARSENED flags if there is new refinement or > coarsening to be done. That would be a much easier library > modification and it would fix the α+A use case, even though it > wouldn't be a step in the right direction long-term. If you do this, > though, you definitely need to add a unit test, so I don't > accidentally break your code again when we ever do get around to long > term refactoring. > --- > Roy -- *~Tim Adowski* |