unit_gradient causing oscillations

[Sean Johnson]

Hey,

I am trying to solve a few coupled hyperbolic equations using discontinuous galerkin finite elements.

In one of my test cases there are oscillations forming. I tracked it down to the gradients of variables that are set to 1 everywhere by:

VectorTools::interpolate(mapping,dof_handler_DG,Functions::ConstantFunction<dim>(1.),alpha_solution);

Obviously changing the final vector for different variables. 

Later, when the gradient of these terms are used the gradient calculated by deal.ii is not 0. I do get very small values on the order of 10^(-14) and when I account for the inverse jacobian factor the unit_gradient returns values ranging between +/- .9x10^(-15) & +/- 2x10^(-15).

These seem very small but due to their values being periodic they quickly grow and destroy stability of my simulation. Are there any known fixes for this?

I haven't noticed any oscillations in the values of the variables from using VectorTools::interpolate for setting initial conditions that would obviously lead to oscillations in the gradient.

Attached is an image on paraview of another variable that depends on one of these oscillating gradients. It has been scaled up by a factor of ~10^(11)

Thanks for your time and any help you can provide,

Sean

[Sean Johnson]

Sorry I realized that image is unhelpful without what it is supposed to be. It should be a straight, flat line of zero values. Now attached is the same image but with the correct value plotted too as a black line.

Best,

Sean

[Wolfgang Bangerth]

Sean:

> VectorTools::interpolate(mapping,dof_handler_DG,Functions::ConstantFunction<dim>(1.),alpha_solution);
>
> Obviously changing the final vector for different variables.
>
> Later, when the gradient of these terms are used the gradient calculated by
> deal.ii is not 0. I do get very small values on the order of 10^(-14) and when
> I account for the inverse jacobian factor the unit_gradient returns values
> ranging between +/- .9x10^(-15) & +/- 2x10^(-15).

Right, and this can not be avoided. Your alpha_solution corresponds to a
function of the form
\sum_j c_j phi_j(x)
where the coefficients c_j in your case happen to all be ones. Then if you
compute the gradient of this function, it is computed as
\sum_j c_j [grad phi_j(x)]
which *should* be zero, but because grad phi_j(x) is not zero, you are adding
together things that are subject to round-off errors and so you have to expect
that the terms in the sum do not exactly cancel but are on the order of
round-off. That's of course exactly what you observe.

> These seem very small but due to their values being periodic they quickly grow
> and destroy stability of my simulation. Are there any known fixes for this?

This is the real problem in your program: That cannot tolerate small
perturbations. This will be true of all perturbations, not just those that
come from round-off. In other words, your algorithm is not *stable*: Small
perturbations grow without bounds. No algorithm that is not stable is useful.
You need to figure out what it is that makes your algorithm instable.

Best
W.


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

[Sean Johnson]

Thanks for the response. I am totally open to the idea of the algorithm not being stable.

The term shown in the graph is from a fairly simple equation. It arises from the need of using the Local Discontinuous Galerkin method for a second order derivative in one of the hyperbolic equations. The equation being solved is:

grad a - q = 0

where a is the variable that is equal to 1 everywhere and q is to be solved for. I try and do a strong form and end up with:

\int_cell  phi * q  =  \int_cell phi * grad a + \int_face phi*( {{a}} - a_interior ) * normal

where {{a}} is the average of the values of a on the faces.

This though leads to the horribly broken solution of q that has positive and negative values and jumps when it should be smooth. Is the answer to use a limiter or filter at this point? I am new to discontinuous galerkin but I haven't seen a filter or limiter being used in local discontinuous galerkin before.

[Wolfgang Bangerth]

On 7/19/24 13:35, Sean Johnson wrote:
>
> The term shown in the graph is from a fairly simple equation. It arises from
> the need of using the Local Discontinuous Galerkin method for a second order
> derivative in one of the hyperbolic equations. The equation being solved is:
> grad a - q = 0
> where a is the variable that is equal to 1 everywhere and q is to be solved
> for. I try and do a strong form and end up with:
>
> \int_cell  phi * q  =  \int_cell phi * grad a + \int_face phi*( {{a}} -
> a_interior ) * normal
>
> where {{a}} is the average of the values of a on the faces.
>
> This though leads to the horribly broken solution of q that has positive and
> negative values and jumps when it should be smooth. Is the answer to use a
> limiter or filter at this point? I am new to discontinuous galerkin but I
> haven't seen a filter or limiter being used in local discontinuous galerkin
> before.

I'm not an expert on DG either, so don't know how to help. But let me ask you
where the the face integral comes from? You're not integrating by parts
anywhere, after all.

[Sean Johnson]

I integrating by parts on the phi * grad a. First time gets it into weak form. Then a second time gets it back into the original equation and into strong form as I understand it. Is this misunderstood?

[Wolfgang Bangerth]

On 7/19/24 14:10, Sean Johnson wrote:
>
> I integrating by parts on the phi * grad a. First time gets it into weak form.
> Then a second time gets it back into the original equation and into strong
> form as I understand it. Is this misunderstood?

That's at least unusual, to integrate twice. I've never seen that done.

[Sean Johnson]

I am following the process from "Nodal Discontinuous Galerkin Methods" by Jan S. Hesthaven and Tim Warburton. It might be particular to discontinuous galerkin to do it this way but I claim to be no expert.

Regardless, if I do the same integral in weak form it still results the same. The equation would change to:

\int_cell  phi * q  =  \int_cell    -grad phi * a + \int_face    phi * {{a}}  * normal

Would an equation so simple need stabilizing?

[Wolfgang Bangerth]

On 7/19/24 14:37, Sean Johnson wrote:
>
> Would an equation so simple need stabilizing?

Yes. It's an advection equation with no diffusion. There needs to be some kind
of stabilization. At the very least you will have to use an upwind flux somewhere.

Best
W>

[Sean Johnson]

Thanks again for your time I really do appreciate it and don't want to waste it. The added equation, although it is an advection equaiont, it is actually representing a diffusion term in the hyperbolic equation. Its just an intermediate term needed for local discontinuous galerkin. Using averages should be stable again based off of the textbook I am using. It is not optimized but it is simpler than doing alternate winding which is the optimized way and to minimize the chance of bugs I am just starting with the simplest form.

I will try and read some papers on stabilizing the method. Thanks!

[Praveen C]

The scheme you wrote looks fine. If a=1, then you should get q=0 or something close to machine zero. What are you getting ? Here you are projecting the gradient onto a discontinuous polynomial space. The blowup probably happens when you solve your advection equation, and maybe you need to look there more carefully.

best

praveen

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/8515c6a2-dfa6-4742-9269-f199200eb1e1n%40googlegroups.com.

[Sean Johnson]

Thanks for the response Praveen,

What I am getting is he image that was in the initial post but here it is again for convenience.

The correct value is in black and is perfectly 0. The calculated values are in white and are on the order of 10^(-12).

The advection equation is solved with a conjugate gradient solver using a multigrid preconditioner as in step 37 of the tutorials since a matrix free method is used.

Best,

Sean Johnson