Advice on using discontinuous methods

General
1

This is a follow up to a previous question I had Seperated Domains within problem . But this question has diverged a bit from the discussion there so I think a seperate post makes sense.

I am trying to solve an advection diffusion problem with a no flux boundary condtion in the middle of the mesh. I was struggling to implement this with mixed dimensional versions of fenics so I am trying to implement it using DG methods. After reading some resources on this method I have managed to get an implementation which is semi-working using a symmetric interior penalty method on the diffusive term, and then balancing the diffusive fluxes against the advective ones.

The solutions appears sensible and the code doesn’t diverge. However, in the system the mass doesn’t appear to be conserved perfectly. In fact if I start off with total mass=4, by time=2 units of simulation time, the mass has gone up to 4.5. If I comment out the advective terms (IE lines 54 and 58-60) then I end up with just simple diffusion with an internal boundary and this seems to work perfectly.

I was wondering if anyone can advise how to improve my simulation with the adevction so that mass is better conserved here? I suspect maybe there is a problem with my formulation of the boundary terms in the symmetric interior penalty method, but I am not sure.

Code:

# Code to solve advection-diffusion  equation using DG method
from fenics import *

set_log_level(30)

T=10
dt=1E-3
theta=0.5

# Create mesh
Lx=1
Ly=1
nx=100
ny=30
mesh = RectangleMesh(Point(-Lx,-Ly),Point(Lx,Ly),nx,ny)

D=0.25
eta=1E-1

# Define function space and functions
V = FunctionSpace(mesh, 'DG', 1)
phi=Function(V)
phi0=Function(V)
v=TestFunction(V)

# initial condition
phi_init=Expression("exp(-pow(x[0]/lx,2)-pow(x[1]/ly,2))",lx=0.2,ly=0.2,degree=2)
# phi_init=Constant(1)
phi0.assign(project(phi_init,V))

# Mark boundary subdomains
class halfway_line( SubDomain ):
  def inside ( self, x, on_boundary ):
    return abs ( x[0] - 0.5 ) < DOLFIN_EPS
  
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
HL=halfway_line()
HL.mark(boundary_parts, 1)

# Set up surface elemtns for fluxes
dS=Measure('dS',domain=mesh,subdomain_data=boundary_parts) # all intenral boundaires
dS_DG=dS(0) # non-physical boundaries due to DG basis
dS_real=dS(1) # internal boundary which is actually real
n=FacetNormal(mesh)  # normal vector, used in fluxes

# Weak form bulk part
phi_star=theta*phi+(1-theta)*phi0 # Crank-Nicolson timestepping
vel=Constant((0.1,0)) # velocity field to advect by
h=mesh.hmin()
F=( 
  (phi-phi0)/dt*v*dx  # time derivative
  + D*dot(grad(phi_star),grad(v))*dx # diffusion
  + div(phi_star*vel)*v*dx # advection
  - D*dot(jump(phi_star,n),avg(grad(v)))*dS_DG # interior flux, this term and next one are from symmetric interior penalty method (SIPG)
  + D*dot(avg(grad(phi_star)),jump(v,n))*dS_DG
  + (eta/avg(h))*dot(jump(phi_star,n),jump(v,n))*dS_DG # interior flux, enforces continuity also from  SIPG
  - phi_star('+')*dot(vel('+'),n('+'))*v('+')*dS_real # balances diffuses flux with advectie one accross internal boundary
  - phi_star('-')*dot(vel('-'),n('-'))*v('-')*dS_real # balances diffuses flux with advectie one accross internal boundary
  - phi_star*dot(vel,n)*v*ds # balances diffuses flux with advectie one accross external boundary
)

# save data
from os.path import expanduser
out_file=expanduser("~/outputs2/DG_test/phi.pvd")
f=File(out_file,"compressed")

for i in range(int(T/dt)):
  
  # Solve nonlinearly as future code will need this
  solve(F==0,phi,
        solver_parameters={"newton_solver":{"relative_tolerance":1E-12,"absolute_tolerance":1E-12,
                                            "maximum_iterations":1000,"relaxation_parameter":1.0}}) # solve using Newton's method

  f<<(phi0,i*dt) # save data
  phi0.assign(phi) # update for next timestep


  total_mass=assemble(phi*dx) # check mass conservation
  print("time: ",i*dt,"total mass: ",total_mass) # print mass conservation
1 Like
2

Are you sure your numerical scheme is consistent and conservative? I wouldn’t expect exact conservation unless it follows from your formulation. If you have a simpler code example, or a reference to a paper or work you’re attempting to replicate, I may have time to check.

3

It’s certainly a consistent and conservative scheme for the diffusion equation. I used the symmetric interior penalty method as described in this playlist 56-58 https://youtube.com/playlist?list=PLMHTjE57oyvpkTPG8ON1w6BChBoapsZTA or in these notes https://www-users.cse.umn.edu/~bcockbur/lecture_notes/DG-4.pdf (page 44 of these notes they discuss these method proposed by Douglas, Jr. and
Dupont in 1976).

However, its not 100% clear to me what modifications are needed when advection also happens, the book says it should be straightforward (page 46 section 5.1) to generalise and so I think I’ve done as they are suggesting, but I’m not fully confident. I should probably read the references in the notes but I was being lazy and hoping somebody here with more experience would see the issues as the notation in the appropriate papers are a bit confusing to me.

If you did have time to look that would be great! I can give a somewhat simpler example (which appears to lose rather than gain mass), but I think some level of complexity is part of the problem itself here to be honest:

# Code to solve advection-diffusion  equation using DG method
from fenics import *

set_log_level(30)

T=10
dt=1E-3
theta=0.5

# Create mesh
Lx=1
Ly=1
nx=100
ny=30
mesh = RectangleMesh(Point(-Lx,-Ly),Point(Lx,Ly),nx,ny)

D=0.25
eta=1E-1

# Define function space and functions
V = FunctionSpace(mesh, 'DG', 1)
phi=Function(V)
phi0=Function(V)
v=TestFunction(V)

# initial condition
phi_init=Expression("exp(-pow(x[0]/lx,2)-pow(x[1]/ly,2))",lx=0.2,ly=0.2,degree=2)
phi0.assign(project(phi_init,V))


# Weak form bulk part
vel=Constant((0.1,0)) # velocity field to advect by
h=mesh.hmin()
n=FacetNormal(mesh)  # normal vector, used in fluxes
F=( 
  (phi-phi0)/dt*v*dx  # time derivative
  + D*dot(grad(phi),grad(v))*dx # diffusion
  + div(phi*vel)*v*dx # advection
  - D*dot(jump(phi,n),avg(grad(v)))*dS # interior flux, this term and next one are from symmetric interior penalty method (SIPG)
  + D*dot(avg(grad(phi)),jump(v,n))*dS
  + (eta/avg(h))*dot(jump(phi,n),jump(v,n))*dS # interior flux, enforces continuity also from  SIPG
  - phi*dot(vel,n)*v*ds # balances diffuses flux with advectie one accross external boundary
)

for i in range(int(T/dt)):
  
  # Solve nonlinearly as future code will need this
  solve(F==0,phi)

  phi0.assign(phi) # update for next timestep


  total_mass=assemble(phi*dx) # check mass conservation
  print("time: ",i*dt,"total mass: ",total_mass) # print mass conservation
4

Why haven’t you considered the inter-element flux arising in the advection term? e.g. Bitbucket. I’m also confused why you’ve left the advective component in the strong form but introduced a facet flux only on the exterior.

Based on the demo you could try the adjustment:

n = FacetNormal(mesh)
un = (dot(vel, n) + abs(dot(vel, n)))/2.0
h=mesh.hmin()
n=FacetNormal(mesh)  # normal vector, used in fluxes
F=( 
  (phi-phi0)/dt*v*dx
  + D*dot(grad(phi),grad(v))*dx
  - inner(phi*vel, grad(v))*dx
  + dot(jump(v), un('+')*phi('+') - un('-')*phi('-') )*dS
  - D*dot(jump(phi,n),avg(grad(v)))*dS
  + D*dot(avg(grad(phi)),jump(v,n))*dS
  + (eta/avg(h))*dot(jump(phi,n),jump(v,n))*dS
  + dot(v, un*phi)*ds
)

But your velocity field of u = (0.1, 0) is inconsistent with the conservation property of the underlying PDE. So you’d need to weakly enforce an outflow condition. In this case you wouldn’t expect to conserve “mass”. If you adjust your velocity field in \Omega = (-1, 1)^2 to something like

x = SpatialCoordinate(mesh)
vel = as_vector((2 * x[1] * (1 - x[0]**2),
                -2 * x[0] * (1 - x[1]**2)))

I’d expect exact conservation to machine precision.

Doing these adjustments and some other changes in your initial condition for aesthetics, I see exact conservation:
abc

1 Like
5

So I’m a little bit confused to be honest. If I go back to adevection-diffusion with CG elements I think the correct weak form is as shown in the following code. I get this by integrating the solely diffusive term by parts, leading to a boundary term of the form D*\nabla( \phi)v \cdot \mathbf{n}ds. But then to generate conservation I need to assume on the boundary that the overall current vanishes on the boundary ie D \nabla \phi=\phi \mathbf{u}. So I replace the boundary term with \phi \mathbf{u} v \cdot \mathbf{n}ds.

This is why in the previous example I keep the advective term in strong from but introduced the flux on the exterior, maybe this is incorrect, but in the case of the CG code shown below the conservation of mass seems to hold up to some numerical error (order 10^-12 at least), which gives me confidence that this works.

This even seems to work with the (0.1,0) velcoity field which is what I would personally expect since for any PDE of the form \frac{ \partial \phi}{ \partial t} + \nabla \cdot J =0, conservation is satisfied if J \cdot \mathbf{n}=0 on the boundary. In this case I am using J=\phi \mathbf{u} - D \nabla \phi, leading to the equality mentioned earlier. This should be valid for all velcoity field \mathbf{u} as far as I know?

This is also the reason why I did not consider the the inter-element flux arising in the advection term. Because I didn’t integrate the advective term by parts, so there is no reason for the boundary term to arise, right? (I think this is most likely where my mistake is, but not sure.)

I’m happy to share my derivation for the weak form I have if it would be helpful to the discussion, its just rather lengthy. Thanks for all your help so far, I can see the weak form you provided is right, I just don’t really understand why mine is wrong.

# Code to solve advection-diffusion  equation using DG method
from fenics import *

set_log_level(30)

T=10
dt=1E-3
theta=0.5

# Create mesh
Lx=1
Ly=1
nx=100
ny=30
mesh = RectangleMesh(Point(-Lx,-Ly),Point(Lx,Ly),nx,ny)

D=0.25
eta=1E-1

# Define function space and functions
V = FunctionSpace(mesh, 'CG', 1)
phi=Function(V)
phi0=Function(V)
v=TestFunction(V)

# initial condition
phi_init=Expression("exp(-pow(x[0]/lx,2)-pow(x[1]/ly,2))",lx=0.2,ly=0.2,degree=2)
# phi_init=Constant(1)
phi0.assign(project(phi_init,V))

# Mark boundary subdomains
class halfway_line( SubDomain ):
  def inside ( self, x, on_boundary ):
    return abs ( x[0] - 0.5 ) < DOLFIN_EPS
  
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
HL=halfway_line()
HL.mark(boundary_parts, 1)

# Set up surface elemtns for fluxes
dS=Measure('dS',domain=mesh,subdomain_data=boundary_parts) # all intenral boundaires
dS_DG=dS(0) # non-physical boundaries due to DG basis
dS_real=dS(1) # internal boundary which is actually real
n=FacetNormal(mesh)  # normal vector, used in fluxes

# Weak form bulk part
phi_star=theta*phi+(1-theta)*phi0 # Crank-Nicolson timestepping
vel=Constant((1,0)) # velocity field to advect by
h=mesh.hmin()
F=( 
  (phi-phi0)/dt*v*dx  # time derivative
  + D*dot(grad(phi_star),grad(v))*dx # diffusion
  + div(phi_star*vel)*v*dx # advection
  - phi_star*dot(vel,n)*v*ds # balances diffuses flux with advectie one accross external boundary
)

# save data
from os.path import expanduser
out_file=expanduser("~/outputs2/DG_test/phi.pvd")
f=File(out_file,"compressed")

for i in range(int(T/dt)):
  
  # Solve nonlinearly as future code will need this
  solve(F==0,phi,
        solver_parameters={"newton_solver":{"relative_tolerance":1E-12,"absolute_tolerance":1E-12,
                                            "maximum_iterations":1000,"relaxation_parameter":1.0}}) # solve using Newton's method

  f<<(phi0,i*dt) # save data
  phi0.assign(phi) # update for next timestep


  total_mass=assemble(phi*dx) # check mass conservation
  print("time: ",i*dt,"total mass: ",total_mass) # print mass conservation