Q&A: Issue with Adaptive Refinement for Higher-Order Elements Q_k+1

[Najwa Alshehri]

Dear developers and users,

I am currently working on a problem involving the coupling of multiple physics. Theoretically, I have shown that a new family of finite element methods is stable for k>= 1. In the case of k = 1, the velocity is approximated by a Q2 finite element [while there are other solutions, I’m focusing here on the velocity for simplicity]. This method should exhibit order "k" convergence,  with adaptive refinement.

The code performs as expected for k=1, both with uniform and adaptive refinement (using an a posteriori error estimator that is theoretically reliable and efficient). However, when moving to k = 2, adaptive refinement does not produce the expected results. I have attached a comparison of the convergence rates between k= 1 (Q2) and k = 2 (Q3). Specifically, while Q2 converges with an H1 norm O(h), I expected Q3 to converge with an H1 norm O(h^2), but the actual result isn’t even proportional to h. Initially, it performs well, but after a few cycles, the issue arises and we lose convergence. (see rate_of_convergance_adp_Q2_Q3.pdf)

Upon further inspection, I found that for k=1 (Q2), the solution (see attached plot of one velocity component) is continuous across the domain ( see Velocity-Q2.png). However, for k = 2 (Q3), discontinuities (jumps) appear at certain hanging nodes ( see Velocity-Q2.png and Velocity-Q3_zoomin.png).

Has anyone encountered this issue before? Could this be related to the handling of hanging node constraints? Even though I suspect that since hanging node I believe built to handle higher order finite elements.

Or is it only an issue with the visualization of higher-order elements? (Note, I did not use GridOut with higher-order mapping in the output.)

I would appreciate any suggestions on how to resolve this issue. 

Best regards,
Najwa

[Najwa Alshehri]

Dear all,

I have an important update on the situation.

I ran step-22 with Q2^d-Q1 and Q3^d-Q2, and I observed the same behavior. Additionally, I found that in my code, even when using Q2, these jumps also appear, though on a smaller scale and in fewer locations.

This brings me to an updated question:

• Is this behavior not significantly affecting the result, and is there any way to improve it?

If this jump is negligible ( or if it is only a visualizing issue), does anyone have any suggestions on where else the issue might lie? The relative errors remain small, but we lose convergence during refinement when using Q3. 

>>> I understand that it is hard to tell without seeing the code, but this question fell into this place after noticing that the jump is not only in my situation.

Notes:

• For refinement, I am using a fixed fraction of 0.3 for refinement and 0.0 for coarsening. I tried different ratios as well.
• I have tried solving using FGMRES with different preconditioners (MG and/or ILU), and the error convergence remains the same, with the jumps still present in the solution plot. The relative error is consistent across different solvers.
• I have correctly applied hanging node constraints at the appropriate stages (DoF setup, after system assembly, and after solving).
• This code should solve for 4 unknowns, 3 of which are vector values, and uses two independent meshes coupled with exact intersection calculation.

If those notes would help anyhow. Thank you all for your support.

Best,
Najwa

[Wolfgang Bangerth]

Najwa:

> I ran step-22 with Q2^d-Q1 and Q3^d-Q2, and I observed the same behavior.
> Additionally, I found that in my code, even when using Q2, these jumps also
> appear, though on a smaller scale and in fewer locations.
>
> This brings me to an updated question:
>

> * Is this behavior not significantly affecting the result, and is there any

> way to improve it?
>
> If this jump is negligible ( or if it is only a visualizing issue), does
> anyone have any suggestions on where else the issue might lie? The relative
> errors remain small, but we lose convergence during refinement when using Q3.

The issue with the gap is likely a visualization artifact:
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#in-my-graphical-output-the-solution-appears-discontinuous-at-hanging-nodes
In other words, what you are seeing is probably not what internally happens.
The documentation of DataOut::build_patches() also discusses this and gives
you a suggestion of how to produce a better visual representation of the
truth. I would try that out first; it is of course also possible that you are
setting up hanging node constraints wrongly, though that's not very likely if
you ask me -- you just have to copy the 6 or 8 lines that deal with hanging
nodes from step-6 and it will work for any scheme.

What it is that affects your convergence rates is of course a different
question, but one I don't think any of us can help with knowing neither the
method nor the implementation in detail :-( Do you get the expected
convergence rate for uniform refinement? That's usually the first place to
start -- make the problem simple -- before using adaptive meshes.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Najwa Alshehri]

Dear Wolfgang,

Thank you for your reply. I 100% agree that this question is not an easy question to ask here in general. The question is meant for the gap that I was think is causing the issue and it shifted to this situation.

Indeed, both cases converge as expected using uniform refinement ( which is quasi-optimal due to the singularity of the solution). Hangging nodes constraints, I believe,  are assigned as I it should be.

It could be an issue related to having one of the meshes circle approximated by a polygon. I have a guess that the error of not fitting the circle exactly is seen worse using higher order elements.  To check this I toke a step back and did the following:

1. I ran another similar problem that is a scaler version of it with no pressure term ( circular mesh is involved) and both Q2 and Q3 converge nicer with adaptive refinement used, however, not optimal for Q3. In this scaler version, we used similar indicators for the adaptive refinement that was also shown to be reliable and efficient. ( Here we did not see bad lose of convergence as before, it is converging as a Q2 which is less than what we expect). ( see plot scalar-circle)
2. I ran a scalar example with no circle involved. So far it seems like it is converging better than a Q2 with adaptive refinement, but still not optimal.( see plot scalar-square). Q2 is converging more than what I expected on the other hand.

Given this, I am trying to run a  new version of the problem of interest with circular mesh is not involved and see how it goes. However, even if this was a factor to our complex issue, I would expect that it become better with refinement which was not the case.

Those are just to update the situation. One more time. I understand that this is not easy to think about with no access to the code. However, since I had the question already posted, I think it would be good to also write what I found ( if I find it  :) ) as the issue with the code. 

In the meantime, any suggestions are highly appreciated.

Best,

Najwa

[Wolfgang Bangerth]

On 9/16/24 05:38, Najwa Alshehri wrote:
>
> It could be an issue related to having one of the meshes circle approximated
> by a polygon. I have a guess that the error of not fitting the circle exactly
> is seen worse using higher order elements.

You may be interested in this discussion:
https://dealii.org/developer/doxygen/deal.II/step_6.html#step_6-Abettermesh

[Najwa Alshehri]

Thanks