### VMS stabilization of 2D steady advection-diffusion fails

151

views

2

Hi,

As I was implementing VMS stabilization of a related problem I ran into some issues. I reduced it to the following MWE. It concerns a simple steady 2D advection-diffusion problem on a square mesh with a diagonally directed advection vector. The problem occurs when I use higher order basis functions (p=2).

The solution completely blows up when the VMS stabilization is included. This can be fixed by switching to SUPG or GLS, but im actually interested in the VMS part. Still, this illustrates that the (sign of the) Laplace operator on the test function is the root of the problem.

- Am I making a mistake with my FEniCS implementation?

- Am I making a mistake with my VMS formulation?

- Is VMS unstuited for the current problem?

Even though I've checked and rechecked the first two points, I can't believe that it would be the third. This MWE should be the simplest example of true VMS (in the sense that p=1 reduces to SUPG and 1D could be reconstructed analytically). I would find it very odd if VMS doesn't work for this case.

I would appreciate any help or advice,

Best,

Stein

As I was implementing VMS stabilization of a related problem I ran into some issues. I reduced it to the following MWE. It concerns a simple steady 2D advection-diffusion problem on a square mesh with a diagonally directed advection vector. The problem occurs when I use higher order basis functions (p=2).

```
from dolfin import *
# Parameters
nu = Constant(0.01) #diffusion coefficient
a_mag = 0.8 #advection magnitude
a = Constant((a_mag*0.5**0.5,a_mag*0.5**0.5)) #advection vector
VMS = False #toggle to include or exclude VMS stabilization
# Create mesh and define function space
mesh = RectangleMesh(Point(-0.5,-0.5),Point(0.5,0.5),10,10)
V = FunctionSpace(mesh, "CG", 2)
# Boundary conditions
def bdy_outer(x,on_boundary):
return on_boundary
bc_outer = Expression("x[0]+x[1]", degree=1,domain=mesh)
bcs = [DirichletBC(V, bc_outer, bdy_outer)]
# Stability parameter
h = CellSize(mesh)
tau = pow( (pow(4.*nu/pow(h,2),2)+pow(2*a_mag/h,2)) ,-0.5)
# Define unknown and test function
v = TestFunction(V)
u = TrialFunction(V)
# Define variational forms
F = -inner(grad(v),a)*u*dx+ nu*inner(grad(u), grad(v))*dx
if VMS:
F += -tau*( -inner(grad(v),a)-nu*div(grad(v)) )*( inner(grad(u),a)-nu*div(grad(u)))*dx
# Solving the problem
u = Function(V)
problem = LinearVariationalProblem(lhs(F),rhs(F), u, bcs)
solver = LinearVariationalSolver(problem)
solver.solve()
# Projecting on a refined mesh for visualization
meshRefine = RectangleMesh(Point(-0.5,-0.5),Point(0.5,0.5),100,100)
V2 = FunctionSpace(meshRefine, "CG", 1)
u2 = project(u,V2)
u2.rename("u", "u")
# Create files for storing results
file = File("solution.pvd")
file << u2
```

The solution completely blows up when the VMS stabilization is included. This can be fixed by switching to SUPG or GLS, but im actually interested in the VMS part. Still, this illustrates that the (sign of the) Laplace operator on the test function is the root of the problem.

- Am I making a mistake with my FEniCS implementation?

- Am I making a mistake with my VMS formulation?

- Is VMS unstuited for the current problem?

Even though I've checked and rechecked the first two points, I can't believe that it would be the third. This MWE should be the simplest example of true VMS (in the sense that p=1 reduces to SUPG and 1D could be reconstructed analytically). I would find it very odd if VMS doesn't work for this case.

I would appreciate any help or advice,

Best,

Stein

Community: FEniCS Project

### 1 Answer

4

This is interesting; it seems like the boundary contributions from the adjoint make a big difference here. (Consider doing the integration-by-parts for \(\mathcal{L}^*\) element-by-element, on the FE space; you'll get extra terms on the element boundaries.) The following gives me nice-looking solutions, even for \(p=5\):

```
# Define variational forms
F = -inner(grad(v),a)*u*dx+ nu*inner(grad(u), grad(v))*dx
if VMS:
#F += -tau*( -inner(grad(v),a)-nu*div(grad(v)) )*( inner(grad(u),a)-nu*div(grad(u)))*dx
def L(u):
return inner(grad(u),a)-nu*div(grad(u))
def Ladj_formal(u):
return -inner(grad(u),a)-nu*div(grad(u))
def badj(u,n):
return nu*inner(grad(u),n)
res = L(u) # - f
uPrime = -tau*res
# Equivalent to what was there originally:
F += Ladj_formal(v)*uPrime*dx
# Boundary terms from deriving the adjoint operator
# in the non-smooth FEM setting:
n = FacetNormal(mesh)
F += badj(v('+'),n('+'))*uPrime('+')*dS \
+ badj(v('-'),n('-'))*uPrime('-')*dS \
+ badj(v,n)*uPrime*ds
```

Hi David,

Thanks for the reply. I like your thought process. Incorporating fine-scale element boundary values into the VMS formulation is part of what I'm working on. I moved away from defining \(u' =\tau \mathcal{R}_{u_h}\) on the boundary though, as I felt that this approximation has little meaning pointwise, and is really only suitable in an integral sense in the element volume. This makes it all the more interesting that the results obtained with your formulation look quite good. One interesting observation in this regard is that the fine-scale is identically zero on the domain boundary (due to the strong BCs), so the 'ds' term should not be there. However, when I remove that term, then the results deteriorate. I'm curious to hear your thoughts on these points.

I came from a spline based implementation when I first observed this problem. In that case, no fine-scale boundary integrals should appear but I still observed the instability. In a coercivity analysis, the negative square Laplacian term is clearly a potential problem. I then find that the stability parameter actually has to be bound from above. Indeed if I multiply \(\tau\) in my initial code by 0.2, then I obtain stable behavior (0.3 being the tipping point, leading to a quite aesthetically pleasing bubbly result). Still, the nonphysical overshoots are not quite mitigated, as is the case in your proposed formulation.

It thus seems that your added terms consistently improve the coercivity of the bilinear form to the point that the original \( \tau \) becomes suitable. This then adds the appropriate amount of artificial diffusion to remove the unphysical boundary layer behavior. It is interesting, though, that coercivity is improved by adding your element boundary terms. I suppose this is largely due to the square gradient terms \( \nabla v\cdot n\, \nabla u \cdot a \), although this can only be a positive contribution on outflow facets. The other \( -\nabla v\cdot n\, \nu \Delta u \) term is also tricky to incorporate in the coercivity analysis.

Thanks for the reply. I like your thought process. Incorporating fine-scale element boundary values into the VMS formulation is part of what I'm working on. I moved away from defining \(u' =\tau \mathcal{R}_{u_h}\) on the boundary though, as I felt that this approximation has little meaning pointwise, and is really only suitable in an integral sense in the element volume. This makes it all the more interesting that the results obtained with your formulation look quite good. One interesting observation in this regard is that the fine-scale is identically zero on the domain boundary (due to the strong BCs), so the 'ds' term should not be there. However, when I remove that term, then the results deteriorate. I'm curious to hear your thoughts on these points.

I came from a spline based implementation when I first observed this problem. In that case, no fine-scale boundary integrals should appear but I still observed the instability. In a coercivity analysis, the negative square Laplacian term is clearly a potential problem. I then find that the stability parameter actually has to be bound from above. Indeed if I multiply \(\tau\) in my initial code by 0.2, then I obtain stable behavior (0.3 being the tipping point, leading to a quite aesthetically pleasing bubbly result). Still, the nonphysical overshoots are not quite mitigated, as is the case in your proposed formulation.

It thus seems that your added terms consistently improve the coercivity of the bilinear form to the point that the original \( \tau \) becomes suitable. This then adds the appropriate amount of artificial diffusion to remove the unphysical boundary layer behavior. It is interesting, though, that coercivity is improved by adding your element boundary terms. I suppose this is largely due to the square gradient terms \( \nabla v\cdot n\, \nabla u \cdot a \), although this can only be a positive contribution on outflow facets. The other \( -\nabla v\cdot n\, \nu \Delta u \) term is also tricky to incorporate in the coercivity analysis.

written
11 weeks ago by
Stein Stoter

3

> I moved away from defining \(u'=\tau\mathcal{R}_{u^h}\) on the boundary though, as I felt that this approximation has little meaning pointwise, and is really only suitable in an integral sense in the element volume.

You might take a look at this old paper by Ken Jansen (if you haven't already):

https://doi.org/10.1016/S0045-7825(98)00284-9

He incorporates the distributional contributions of the residual on interior facets.

> One interesting observation in this regard is that the fine-scale is identically zero on the domain boundary (due to the strong BCs), so the 'ds' term should not be there. However, when I remove that term, then the results deteriorate. I'm curious to hear your thoughts on these points.

While it is true that the exact \(u'\) must be zero to satisfy \(\overline{u} + u' = 0\) on the boundary when strongly enforcing \(\overline{u}=0\) there, the substitution \(u'\gets-\tau\mathcal{R}\) can be nonzero, so, in the discrete setting, there is a difference.

You might take a look at this old paper by Ken Jansen (if you haven't already):

https://doi.org/10.1016/S0045-7825(98)00284-9

He incorporates the distributional contributions of the residual on interior facets.

> One interesting observation in this regard is that the fine-scale is identically zero on the domain boundary (due to the strong BCs), so the 'ds' term should not be there. However, when I remove that term, then the results deteriorate. I'm curious to hear your thoughts on these points.

While it is true that the exact \(u'\) must be zero to satisfy \(\overline{u} + u' = 0\) on the boundary when strongly enforcing \(\overline{u}=0\) there, the substitution \(u'\gets-\tau\mathcal{R}\) can be nonzero, so, in the discrete setting, there is a difference.

written
11 weeks ago by
David Kamensky

Please login to add an answer/comment or follow this question.

https://doi.org/10.1016/S0045-7825(98)00079-6