Major error with AMR mesh for DG problems

[class4kayaker]

When a mesh boundary is refined using AMR, the ties to periodic neighbors appear to be lost. Note that the problem does not appear to occur if the refinement is done during the pre-refinement step, but does occur when refined later.

I would need to look closer at the problem/look deeper into the implementation to identify the source of the issue, but the nature of the issue results in significantly incorrect behavior.

The attached file is a test which advects a simple composition band band on a periodic domain in a manner which can be expected to cause changes in boundary refinements and demonstrates the issue.

test_periodic.prm.txt

[yinghe616]

@class4kayaker During hackathon, I only checked the DG periodic case with a uniform mesh but not on ARM, and also only checked AMR algorithm with composition approximate gradient criterion when using discontinuous composition (doesn't have to be used for DG method) for non-periodic cases. Yes, it seems the AMR mesh doesn't work well for DG method based on the parameter file you provided. I think the possible reason could be that after the mesh refined, during resetting the dof process, current code only updated the constraint matrix on the periodic boundaries used for FEM, however DG method doesn't use constraint matrix to enforce the periodic boundaries. I need to look into core.cc for solutions, but thanks for pointing this out. @bangerth any comments?

[yinghe616]

@class4kayaker I looked at the code, it seems to me that the periodic boundary setting is fine for the DG method; then I tried to use finer mesh and a larger value refinement fraction number, I can see some nonzero values appear at the periodic boundary in the top of the box at a later time. I think it verified that the codes at least can recognize the periodic boundary and assign some values. However I found that the values in the top are quire large, it seems a large overshoot. Then I used the DG limiter, it works well now.

Some changes I made on the prm file:

subsection Mesh refinement set Minimum refinement level = 3 set Initial adaptive refinement = 1 set Initial global refinement = 5 set Time steps between mesh refinement = 1 set Coarsening fraction = 0.05 set Refinement fraction = 0.95 set Strategy = composition approximate gradient end

subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters set Use limiter for discontinuous composition solution = true # apply the limiter to the DG solutions set Global composition maximum = 1.0 set Global composition minimum = 0.0 end end

[bangerth]

@class4kayaker -- I don't have the time right now to look at it myself, but as a starting point for anyone who wants to look into this: can you attach pictures to this issue that show what is happening? Oftentimes, looking at pictures already gives a hint of what may be happening, and we may be able to suggest a place to go look next.

[class4kayaker]

Sure, the issue is that for whatever reason a periodic boundary is not being used after the boundary has mesh adjustments despite being specified. The attached parameter file does not produce the relevant error using a continuous composition field, but if changed to use a discontinuous composition field the error occurs. I noted some additional suggestive pieces of information based on the problem I first saw the issue occur (a periodic AMR VoF test problem similar the the problem I included a parameter file for).

Ying, I do not believe the issue is with your code specifically, but rather with the AMR mesh update somehow not preserving the periodicity data. I originally noted the problem when testing my VoF code for periodic boundary conditions and attempted to identify the issue with that refinement code, but did not think to test the same type of problem with a DG composition condition until today. The fact that both algorithms had the same issue suggests that the problem is not with the composition code, but rather with the mesh data. The images should highlight what problem I am referring to.

The test problem is a single C=1 band being advected in the negative-Y direction using a box geometry with periodic-Y boundary conditions. As you can see, when the composition field actually reaches the lower boundary, the boundary does not appear to be treated as periodic. When testing on a VoF periodic AMR problem (how I originally identified the issue), it appeared that the bug only occurred if the boundary in question had the refinement level change during the test run, in particular if the mesh was refined on that boundary.

Initial condition (t=0)

Image during the problem (t=1/2 period)

Final configuration (t= full period)

The periodic boundary being treated as an outflow boundary condition is obvious from the above screenshots.

My current guess at the source of the error is that the periodic neighbor specifications are lost when refining a periodic boundary conditions. I am not familiar enough with the AMR code to guess where to start looking.

[bangerth]

Thanks for the picture. It's clear something is wrong.

The periodicity in deal.II is handled at the level of the triangulation. If there was a bug there, it would also show up if you had a continuous finite element. Did you try that? Otherwise, it must have something to do with the construction of DG fluxes.

[class4kayaker]

As I noted, I tested with a continuous finite element, and the problem did not occur. I have attached second test problem using standard Aspect which will show the same behavior from parts of the VoF test run that suggested the issue is not with the flux construction code (assembly) but likely to be an issue with the triangulation data. In particular, when the boundary is initially refined and unchanged the algorithm does not have any issues (first period), but upon a change in the boundary (second period) the periodic boundary condition fails to behave correctly.

test_periodic.prm.txt

(uses Ying's suggested modifications for better behavior)

Initial conditions (upper layer necessary to maintain high refinement of boundary, could be part of issue)

Almost 1 period of time

During second advection over periodic boundary

Final state

It might be the flux calculation, but if so the most likely source is an error/undocumented caveat to periodic neighbor methods.

[yinghe616]

For continuous element, the periodicity is enforced strongly by the constrained matrix (so the inconsistent of refinement levels between top and bottom boundaries may not cause any trouble). For discontinuous element, the periodicity is enforced weakly by flux calculations, upwind for this case (so the inconsistent refinement levels may cause such behavior. Note, the cells near the top boundary should be refined when the banded box touches the lower boundary, which is the direction the banded box is going to move.) I think to further debug this issue, one can try to enforce the cell refinement levels to be the same for both the top and bottom periodic boundaries. On Thu, Jun 22, 2017 at 8:36 PM Jonathan Robey notifications@github.com wrote:

As I noted, I tested with a continuous finite element, and the problem did not occur. I have attached second test problem using standard Aspect which will show the same behavior from parts of the VoF test run that suggested the issue is not with the flux construction code (assembly) but likely to be an issue with the triangulation data. In particular, when the boundary is initially refined and unchanged the algorithm does not have any issues (first period), but upon a change in the boundary (second period) the periodic boundary condition fails to behave correctly.

test_periodic.prm.txt https://github.com/geodynamics/aspect/files/1096693/test_periodic.prm.txt

(uses Ying's suggested modifications for better behavior)

Initial conditions (upper layer necessary to maintain high refinement of boundary, could be part of issue) [image: advect-init] https://user-images.githubusercontent.com/1050493/27465615-820429f4-5789-11e7-949a-65b7db26cb90.png

Almost 1 period of time [image: advect-1] https://user-images.githubusercontent.com/1050493/27465612-8202be0c-5789-11e7-8458-e6b710e9e5b9.png

During second advection over periodic boundary [image: advect-problem] https://user-images.githubusercontent.com/1050493/27465613-8202e436-5789-11e7-8bda-6ce8839de976.png

Final state [image: advect-final] https://user-images.githubusercontent.com/1050493/27465614-8203d7f6-5789-11e7-88d9-50d6d79a4c66.png

It might be the flux calculation, but if so the most likely source is an error/undocumented caveat to periodic neighbor methods.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/geodynamics/aspect/issues/1822#issuecomment-310562522, or mute the thread https://github.com/notifications/unsubscribe-auth/AL9LHNCeIjc9-I7oL9yCJztCnWf1mM7Nks5sGzKlgaJpZM4OC8Dj .

[yinghe616]

@class4kayaker @bangerth there is indeed a bug in the DG flux construction on the AMR grid. See detail in the PR #1823 I'll continue to work on fixing it. Thanks @class4kayaker to bring this question.

[class4kayaker]

Looking into the VoF implementation, and it looks like Ying has correctly diagnosed the issue. I believe I was somewhat confused due to the initial refinement sidestepping the issue by maintaining consistent refinement levels.

[bangerth]

This relates to #1522.

[bangerth]

@class4kayaker -- just as a follow-up to explain why I knew that the problem must be in how we compute the flux terms: If the problem does not exist with continuous elements, then apparently the triangulation got it right in marking which cells are periodic neighbors of which other cells. Because if the triangulation had gotten that part wrong, then the continuous element case could not possibly have worked because it would have worked on incorrect data. So the problem can't have been in the triangulation.

If, however, the triangulation is correct in stating which cells are periodic neighbors to which other cells, then the problem must be somewhere in the difference between how we handle continuous and discontinuous elements. The former use constraints, which work, whereas the latter uses flux terms, which don't work. My conclusion is that the problem must be with the flux terms -- as indeed @yinghe616 quickly demonstrated.

I don't know if there is a lesson here that I can clearly articulate, but the point is that oftentimes asking the right questions can get you pretty quickly to identifying where the problem may be. It typically comes down to figuring out whether "I know that this condition ought to be true" is really correct -- in the current case, the condition is "if the periodicity information is wrong, then the continuous element cannot work correctly". The second half of the statement is easy to test (as you did), and because it is wrong, the assumption in the first half of the statement cannot have been correct.

[gassmoeller]

Closed by #2132