FEniCS Project

DG upwind crashes in 2D

Asked by Michał Wichrowski

I have working 1D code, it works fine:
-------------------------------------------
from dolfin import *

# Load mesh
mesh = UnitInterval(14)
# Defining the function spaces
V_dg = FunctionSpace(mesh, "DG", 8)
V_cg = FunctionSpace(mesh, "CG", 5)

V_u = FunctionSpace(mesh, "CG", 2)

# Create velocity Function

vel= Expression('(x[0]-0.5)*2')

u = project(vel, V_u)
plot(u)

# Test and trial functions
v = TestFunction(V_dg)
phi = TrialFunction(V_dg)

# Diffusivity
kappa = Constant(8E-8)

# Source term
f1 = Expression("pi*2*(x[0]-0.5)*sin(pi*2*(x[0]-0.5))")
f2 = kappa*pi*pi*Expression("cos(pi*2*(x[0]-0.5) )")
f = f1+f2
F = project(f, V_cg)
plot(F)

# Penalty term
alpha = Constant(5.0)
gamma = Constant(12.0)

# Mesh-related functions
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2

fi = -(dot(u, n) - abs(dot(u, n)))/(2.0*dot(u,n))

boundary = Expression('-2+2*x[0]')

sigma = - div(u)
#sigma = Constant('0')

# Bilinear form
a_int = dot(grad(v), kappa*grad(phi))*dx + dot(grad(v), u)*phi*dx - v*phi*sigma*dx

a_fac = 0*kappa('+')*(alpha('+')/h('+'))* dot(jump(v, n), jump(phi, n))*dS \
      - dot(avg(kappa*grad(v)), jump(phi, n))*dS \
      + dot(jump(v, n), avg(kappa*grad(phi)))*dS

a_vel = dot((v('+')*n('+') + v('-')*n('-')),(u('+')*phi('+')*fi('+') + u('-')*phi('-')*fi('-')))*dS

a_bound = gamma/h*v*phi*ds
#a_bound = dot(kappa*grad(v),n)*phi*ds - v*dot(kappa*grad(phi),n)*ds

a = a_int + a_fac + a_vel + a_bound

# Linear form
L = v*f*dx + gamma/h*v*boundary*ds
#L = v*f*dx + dot(kappa*grad(v),n)*boundary*ds - v*boundary*dot(u,n)*fi*ds

# Solution function
phi_h = Function(V_dg)

# Assemble and apply boundary conditions
A = assemble(a)
b = assemble(L)

# Solve system
solve(A, phi_h.vector(), b)

# Plot solution
plot(phi_h)
interactive()
-------------------------------
I've changed mesh and velocity:

mesh = UnitSquare(30,60)
V_u = VectorFunctionSpace(mesh, "CG", 2)
vel = Expression (('x[0]-0.5','0'))
 and I've got error massege:
---------------------
    plot(phi_h)
  File "/usr/lib/python2.7/dist-packages/dolfin/common/plotting.py", line 121, in dolfin_plot
    return viper_dolfin.plot(make_viper_object(object, mesh=mesh), *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/viper/viper_dolfin.py", line 452, in plot
    fig = _plotter.plot(data, **kwargs)
  File "/usr/lib/python2.7/dist-packages/viper/viper_dolfin.py", line 401, in plot
    return self.autoplot(plot_object, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/viper/viper_dolfin.py", line 415, in autoplot
    plotter = Viper(plot_object, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/viper/viper_dolfin.py", line 60, in __init__
    self.plot(data, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/viper/viper_dolfin.py", line 85, in plot
    plot_method(data, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/viper/viper_dolfin.py", line 153, in plot_genericfunction
    assert vmax >= vmin, "Empty range, please specify vmin and/or vmax"
AssertionError: Empty range, please specify vmin and/or vmax
----------------------

Question information

Language:
English Edit question
Status:
Solved
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Solved by:
Michał Wichrowski
Solved:
Last query:
Last reply:
Revision history for this message
Michał Wichrowski (mtwichrowski) said : 		#1

I've solved the problem:
fi = -(dot(u, n) - abs(dot(u, n)))/(2.0*dot(u,n))
should be changed to:
fi = conditional (ge(dot(u,n),0),1,0)

Revision history for this message
Derek Monahan (derek-d-monahan) said : 		#2

Hi Michal,

As I am currently still learning FEniCS I am trying to reverse engineer your posted code to get an understading of DG implementation. It seems like an ideal learning example as it's non-trivial but short. Your solution looks a bit like the solution to the "Hemker Problem" (Donea & Huerta, "Finite Element Methods for Flow Problems", p78) but not quite the same. Could I ask if you could specifiy the exact problem you are solving?

Thanks,

Derek

Revision history for this message
Michał Wichrowski (mtwichrowski) said : 		#3

Hi Derek
It was Hemker problem in 1D, I've scaled x variable to fit unit
interval. That's why 2*(x[0]-0.5) appears everywhere instead of x[0].
I've made a mistake somewhere, so that the solution is little
different from solution of Hemker problem. Probably in source therm f1
should be multiplied by 2.

I was testing DG scheme for advection-diffusion problems with
discontinuities. Hemker problem was a good test case in 1D but I've
needed something like that in 2D and 3D. The simplest idea was to use
the same equation so code could used with very little changes. The
line of discontinuity could be easily changed by changing velocity, so
the program could be used to demonstrate what happens if discontinuity
is not on face. This problem can also demonstrate how important is to
use conservative scheme - when I've used method from
undocumented/dg-advection-diffusion (in 1D) i've get 2
discontinuities - at x=0 and x=1.

 I've made a mistake somewhere, so that the solution is little
different from solution of Hemker problem. Probably in source therm f1
should be multiplied by 2.
Best,
Michał

Revision history for this message
Derek Monahan (derek-d-monahan) said : 		#4

Thanks for the detailed and quick reply Michal,

I've tried your correction and that does indeed seem to give much better agreement. This is a great example of he DG method and superbly demonstrates it's virtues over stabilization schemes such as SUPG. If you have a working 2d example I would be very grateful if you could post it, as I said I have tried your suggested changes but I get an all zero solution.

I'll be playing with your code for a few hours so I'll let you know if I find anything interesting :)

Thanks again,

Derek