Slope limiting for DG

variational formulation
1

Hi,

has anyone attempted to use slope limiters for DG schemes in fenics? Right now I’m solving a 2D euler-like system.

Is there any good references on slope limiters for DG methods or projects that use them?

There’s a description of a slope limiting procedure in this reference but right now I wouldn’t know how to start coding it.

Thanks in advance.

2

In my experience, “good” slope limiters are a real pain to implement in dolfin since you need to increase the compact support of bases to extend further than the local neighbours of cells.

Perhaps consider the work in Ocellaris, e.g. here.

1 Like
3

Thanks, I will check it.

At first glance the code looks quite complex indeed.

Do you know of any other options to tame problems arising from negative density solutions due to oscillations in higher order discretisations?

4

The (little) work I’ve done on approximations of non-smooth compressible flow solutions was not with FEniCS. I got around it with adaptive hp-refinement as developed by my PhD advisor. I’d resolve regions with large transitions or shocks with low order approximations.

With the new developments in basix I really hope hp-refinement could be in dolfinx’s future. But this requires a research question to make it worth the time to implement.

You could also consider artificial viscosity, e.g. eq (14) of here.

1 Like
5

Great, thanks for all the information.

6

Hi, @nate, I’m having some problems trying to implement the artificial viscosity solution. Namely, when I calculate any differential operator on the flux tensor I get a Zero object, for example:

import dolfin as df
import ufl as ufl
import numpy as np

mesh = df.UnitSquareMesh(32, 32)

V = df.VectorFunctionSpace(mesh, 'DG', 0, 3)

u = df.Expression(('sin(pi*x[0]) * sin(pi*x[1])',
                        '0.',
                        '0.'), pi = np.pi, degree=2)

u = df.project(u, V)

def F(u):
    return df.as_tensor([[u[1], u[2]],
                         [u[1]*u[1]/u[0], u[1]*u[2]/u[0]],
                         [u[1]*u[2]/u[0], u[2]*u[2]/u[0]]])

F = F(u)

div_F = ufl.div(F)

print(div_F)

outputs:

0 (shape (3,))

This doesn’t happen when i calculate derivatives of test functions in the convective part of the flux:

v = df.TestFunction(V)
print(ufl.grad(v))

Output:

grad(v_0)

Am I doing something wrong in the calculation of the diffusive flux?

7

Looks like your function space V_\mathrm{DG}^{h,p} is degree p=0, so \mathrm{div}(F) = 0 \quad \forall F \in V_\mathrm{DG}^{h,0}. A first order accurate DG method should be highly dissipative as-is and not require stabilisation.

1 Like
8

True, my bad. I was testing everything with the order zero discretization which indeed works well and I did not realize the gradients vanish there. Thanks for the quick reply.

9

Hi @nate, I just remembered this topic and wanted to let you know that I could make the artificial viscosity solution work and it does work quite well so thank you for the indications.

I post a couple (academic) examples in case anybody is interested:

  1. The first is the stabilization of an isentropic shock tube with order 1 polynomials when introducing artificial viscosity.
  2. The second is a supersonic (mach 3) forward facing step on a very coarse mesh.

shock_tube_stabilization
shock_tube_stabilization1066×461 34.5 KB

forward_facing_step.PNG
forward_facing_step.PNG891×310 108 KB

1 Like
10

These look amazing, thanks for sharing. Is your code open or your work in a preprint or paper somewhere?

1 Like
11

I would like to make it open at some point but it will take some time to clean and document it. I’ll keep the forum informed anyway.

2 Likes