How can I write a simple upwind scheme using DG0 and RT1 elements?

Question

I am trying to construct a simple upwind scheme, using DG0 and RT1 elements.

What is the best way to do this? I have managed to get something that seems to work,
but it seems a very inelegant way to perform a simple calculation.

from dolfin import *                         

# mesh with one internal boundary                                                                                                        
mesh = UnitSquareMesh(1,1)                                                                              
RT = FunctionSpace(mesh, "RT", 1)                                                                       
DG = FunctionSpace(mesh, "DG", 0)                                                                       

# velocity                                                                                                        
U = Function(RT)                                                                                        

# advected quantity (1 in left cell, 2 in right cell)                                                                                         
C = Function(DG)                                                                                        
C.interpolate(Expression("(x[0]>0.5)+1.0"))                                                             

n = FacetNormal(mesh)                                                                                   
un = dot(U,n)                                                                                           
up = (un + abs(un))/2.0                                                                                 
um = (un - abs(un))/2.0                                                                                 

Ct = TestFunction(DG)                                                                                   
flux = ((up*Ct*C)('+') + (up*Ct*C)('-') + (um*Ct)('+')*C('-') + (um*Ct)('-')*C('+'))*dS                 

U.interpolate(Constant((0.1,0.0)))                                                                      
print assemble(flux).array()                                                                            

U.interpolate(Constant((-0.1,0.0)))                                                                     
print assemble(flux).array()  

Answer 1

I think this is an equivalent expression, which is simpler:

flux = dot(jump(Ct), jump(up*C))

where C is a DG0 Function with the field to be advected, Ct is a TestFunction
and "up" is the upwind velocity as defined above...

This is more-or-less what is in the dg-advection-diffusion demo...

Comments

What is C0?

---

Updated answer...

Answer 2

Perhaps you can use the same trick as in the dg-advection-diffusion demo:

un = (dot(u, n) + abs(dot(u, n))) / 2

Comments

Yes, that's how I started out.

The outward flux can be calculated as:

$$ Q_{out} = (u\cdot \hat{n} + || u\cdot \hat{n}||) / 2 $$

but when that is restricted on the facet, it (obviously) only gives a non-zero value in the cell with the outflow. That is fine for updating the cell with the outflow, but then you have to consider the inflow separately. I was just wondering if anyone had anything simpler than what I had written above (which does work).

Also, the restriction operators are a bit 'restrictive' - it might be useful to define ('in') or ('interior') and ('out') or ('exterior') for the Facets which are inside or outside the cell, also ('all') for all Facets, and maybe an alternation operator, e.g.

alt(A, B) = A('+')*B('-') + A('-')*B('+')

This illustrates the simple behaviour of outflow:-

n = FacetNormal(mesh)
un = dot(U,n)
up = (un + abs(un))/2.0

Ct = TestFunction(DG)

print "u = (0.1, 0.0)"
U.interpolate(Constant((0.1, 0.0)))
print assemble((up*Ct)('-')*dS).array()
print assemble((up*Ct)('+')*dS).array()
print
print "u = (-0.1, 0.0)"
U.interpolate(Constant((-0.1, 0.0)))
print assemble((up*Ct)('-')*dS).array()
print assemble((up*Ct)('+')*dS).array()

outputs:

u = (0.1, 0.0)
[ 0.   0.1]
[ 0.  0.]

u = (-0.1, 0.0)
[ 0.  0.]
[ 0.1  0. ]

---

Hi Logg, what is the reference paper for dg-advection-diffusion demo?

---

Hi. everyone:
I think it may be helpful for me to paste my problem here. I am trying to write a simple case for advection of Guass Function with DG method with inlet and outlet. I tried on a unitsquare mesh and try to mimic what in the dg-advection-diffussion demo in constructing flux as you guys discuss before, where I define the adevction velocity on a CG space.I am now trying Eular time stepping,which gives diverged after several steps.here is the code:

from dolfin import *
mesh = UnitSquare(32,1)
c=-1
c0=Expression('exp(-pow((x[0]+a*t),2)/0.025)',a=c,t=0)
def c0_boundary(x):
    tol = 1E-15
    return abs(x[0]) < tol


# Create function space constrained by periodic BC
V_dg = FunctionSpace(mesh, 'DG',1)
V_cg = VectorFunctionSpace(mesh,'CG',1)

#define the advection velocity
u0=Expression(('a','0'),a=c)
U=interpolate(u0,V_cg)

# Define normal component, mesh size and right-hand side
n = FacetNormal(mesh)
un = dot(U,n)
up = (un + abs(un))/2.0

c=TrialFunction(V_dg)
v=TestFunction(V_dg)

a_M = c*v*dx
M=assemble(a_M)

c=Function(V_dg)

k=[]
for i in range(4):
    k.append(Function(V_dg))

T = 0.2
t = 0
dt =0.001
alpha=0
while t<=T:
    a_int=-dot(grad(v),U*c)*dx 
    a_flux=dot(jump(v), un('+')*c('+') - un('-')*c('-') )*dS  + dot(v, un*c)*ds
    a=a_int + a_flux
    b=assemble(a)

    c0.t=t
    bc=DirichletBC(V_dg,c0,c0_boundary)
    bc.apply(M,b)
    solve(M,k[0].vector(),b)

    c.vector()[:]=k[0].vector()
    t += dt
    plot(c,rescale=True)
 can anyone  tell me what is wrong here