Re: [Libmesh-devel] EquationSystems reinit() without refine/coarsen calls

[Roy S. <roy...@ic...>]

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