Stabilized convection-difusion

Convection diffusion equation

Find approximate solution to following linear PDE

\[\begin{split}u_t + \mathbf{b}\cdot\nabla{u} - \operatorname{div}(K \nabla u) &= f \quad\text{ in }\Omega\times[0, T], \\ u &= u_\mathrm{D} \quad\text{ in }\Gamma\times[0, T], \\ u &= u_0 \quad\text{ on }\Omega\times{0} \\\end{split}\]

using \(\theta\)-scheme discretization in time and arbitrary FE discretization in space with data

• \(\Omega = [0, 1]^2\)
• \(T = 10\)
• \(K = \frac{1}{\mathrm{Pe}}\)
• \(\mathbf{b} = \left( -(y-\tfrac{1}{2}), x-\tfrac{1}{2}\right)\)
• \(f = \vec{0}\)
• \(u_0(\mathbf{x}) = \left( 1 - 25 \operatorname{dist}\left(\mathbf{x}, \left[\frac{1}{4}, \frac{1}{4}\right]\right) \right) \chi_{ B_{1/5}\left(\left[\frac{1}{4}, \frac{1}{4}\right]\right) }\)
• \(u_\mathrm{D} = 0\)

where \(\chi_X\) is a characteristic function of set \(X\), \(B_R(\mathbf{z})\) is a ball of radius \(R\) and center \(\mathbf{z}\) and \(\operatorname{dist}(\mathbf{p},\mathbf{q})\) is Euclidian distance between points \(\mathbf{p}\), \(\mathbf{q}\).

Reference solution

• simple straight-forward Galerkin FEM

  import time
  import os
  import math
  from dolfin import *
  
  # get file name
  fileName = os.path.splitext(__file__)[0]
  
  parameters['form_compiler']['cpp_optimize'] = True
  parameters['form_compiler']['optimize'] = True
  
  # Parameters
  Pe = Constant(1e10)
  t_end = 10 
  dt = 0.1
  
  # Create mesh and define function space
  mesh = RectangleMesh(0, 0, 1, 1, 40, 40, 'crossed')
  
  # Define function spaces
  V = FunctionSpace(mesh, "CG", 1)
  
  #ic= Expression("((pow(x[0]-0.25,2)+pow(x[1]-0.25,2))<0.2*0.2)?(-25*((pow(x[0]-0.25,2)+pow(x[1]-0.25,2))-0.2*0.2)):(0.0)")
  ic= Expression("((pow(x[0]-0.3,2)+pow(x[1]-0.3,2))<0.2*0.2)?(1.0):(0.0)", domain=mesh)
  
  b = Expression(("-(x[1]-0.5)","(x[0]-0.5)"), domain=mesh)
  
  bc=DirichletBC(V,Constant(0.0),DomainBoundary())
  	
  # Define unknown and test function(s)
  v = TestFunction(V)
  u = TrialFunction(V)
  
  u0 = Function(V)
  u0 = interpolate(ic,V )
  
  # STABILIZATION
  h = CellSize(mesh)
  n = FacetNormal(mesh)
  theta = Constant(1.0)
  
  # Define variational forms
  a0=(1.0/Pe)*inner(grad(u0), grad(v))*dx + inner(b,grad(u0))* v *dx
  a1=(1.0/Pe)*inner(grad(u), grad(v))*dx + inner(b,grad(u))* v *dx
  
  A = (1/dt)*inner(u, v)*dx - (1/dt)*inner(u0,v)*dx + theta*a1 + (1-theta)*a0
  
  F = A
  
  # Create files for storing results
  file = File("results_%s/u.xdmf" % (fileName))
  
  u = Function(V)
  ffc_options = {"optimize": True, "quadrature_degree": 8}
  problem = LinearVariationalProblem(lhs(F),rhs(F), u, [bc], form_compiler_parameters=ffc_options)
  solver  = LinearVariationalSolver(problem)
  
  
  u.assign(u0)
  u.rename("u", "u")
  
  # Time-stepping
  t = 0.0
  
  file << u
  
  while t < t_end:
  
  	print "t =", t, "end t=", t_end
  
  	# Compute
  	solver.solve()
          plot(u)
  	# Save to file
  	file << u
  
  	# Move to next time step
  	u0.assign(u)
  	t += dt
  

• with SUPG stabilization

  import time
  import os
  import math
  from dolfin import *
  
  # get file name
  fileName = os.path.splitext(__file__)[0]
  
  parameters['form_compiler']['cpp_optimize'] = True
  parameters['form_compiler']['optimize'] = True
  
  # Parameters
  Pe = Constant(1e10)
  t_end = 10 
  dt = 0.1
  
  # Create mesh and define function space
  mesh = RectangleMesh(0, 0, 1, 1, 40, 40, 'crossed')
  
  # Define function spaces
  V = FunctionSpace(mesh, "CG", 1)
  
  #ic= Expression("((pow(x[0]-0.25,2)+pow(x[1]-0.25,2))<0.2*0.2)?(-25*((pow(x[0]-0.25,2)+pow(x[1]-0.25,2))-0.2*0.2)):(0.0)")
  ic= Expression("((pow(x[0]-0.3,2)+pow(x[1]-0.3,2))<0.2*0.2)?(1.0):(0.0)", domain=mesh)
  
  b = Expression(("-(x[1]-0.5)","(x[0]-0.5)"), domain=mesh)
  
  bc=DirichletBC(V,Constant(0.0),DomainBoundary())
  	
  # Define unknown and test function(s)
  v = TestFunction(V)
  u = TrialFunction(V)
  
  u0 = Function(V)
  u0 = interpolate(ic,V )
  
  # STABILIZATION
  h = CellSize(mesh)
  n = FacetNormal(mesh)
  theta = Constant(1.0)
  
  nb = sqrt(inner(b,b))
  tau = 0.5*h*pow(4.0/(Pe*h)+2.0*nb,-1.0)
  
  # first alternative: redefine the test function
  #v=v+tau*inner(b,grad(v))
  
  # second alternative: explicitly write the additional terms
  r=((1/dt)*(u-u0) + theta*((1.0/Pe)*div(grad(u)) + inner(b,grad(u)))+ (1-theta)*((1.0/Pe)*div(grad(u0)) + inner(b,grad(u0))) )* tau * inner(b,grad(v))*dx
  
  # Define variational forms
  a0=(1.0/Pe)*inner(grad(u0), grad(v))*dx + inner(b,grad(u0))* v *dx
  a1=(1.0/Pe)*inner(grad(u), grad(v))*dx + inner(b,grad(u))* v *dx
  
  A = (1/dt)*inner(u, v)*dx - (1/dt)*inner(u0,v)*dx + theta*a1 + (1-theta)*a0
  
  F = A + r
  
  # Create files for storing results
  file = File("results_%s/u.xdmf" % (fileName))
  
  u = Function(V)
  ffc_options = {"optimize": True, "quadrature_degree": 8}
  problem = LinearVariationalProblem(lhs(F),rhs(F), u, [bc], form_compiler_parameters=ffc_options)
  solver  = LinearVariationalSolver(problem)
  
  
  u.assign(u0)
  u.rename("u", "u")
  
  # Time-stepping
  t = 0.0
  
  file << u
  
  while t < t_end:
  
  	print "t =", t, "end t=", t_end
  
  	# Compute
  	solver.solve()
          plot(u)
  	# Save to file
  	file << u
  
  	# Move to next time step
  	u0.assign(u)
  	t += dt
  

• with IP (interior penalty) approximation

  import time
  import os
  import math
  from dolfin import *
  
  # get file name
  fileName = os.path.splitext(__file__)[0]
  
  parameters['form_compiler']['cpp_optimize'] = True
  parameters['form_compiler']['optimize'] = True
  
  # Parameters
  Pe = Constant(1e10)
  t_end = 10 
  dt = 0.1
  
  # Create mesh and define function space
  mesh = RectangleMesh(0, 0, 1, 1, 40, 40, 'crossed')
  
  # Define function spaces
  V = FunctionSpace(mesh, "CG", 1)
  
  #ic= Expression("((pow(x[0]-0.25,2)+pow(x[1]-0.25,2))<0.2*0.2)?(-25*((pow(x[0]-0.25,2)+pow(x[1]-0.25,2))-0.2*0.2)):(0.0)")
  ic= Expression("((pow(x[0]-0.3,2)+pow(x[1]-0.3,2))<0.2*0.2)?(1.0):(0.0)", domain=mesh)
  
  b = Expression(("-(x[1]-0.5)","(x[0]-0.5)"), domain=mesh)
  
  bc=DirichletBC(V,Constant(0.0),DomainBoundary())
  	
  # Define unknown and test function(s)
  v = TestFunction(V)
  u = TrialFunction(V)
  
  u0 = Function(V)
  u0 = interpolate(ic,V )
  
  # STABILIZATION
  h = CellSize(mesh)
  n = FacetNormal(mesh)
  theta = Constant(1.0)
  
  alpha = Constant(0.1)
  
  # Define variational forms
  a0=(1.0/Pe)*inner(grad(u0), grad(v))*dx + inner(b,grad(u0))* v *dx
  a1=(1.0/Pe)*inner(grad(u), grad(v))*dx + inner(b,grad(u))* v *dx
  
  A = (1/dt)*inner(u, v)*dx - (1/dt)*inner(u0,v)*dx + theta*a1 + (1-theta)*a0
  
  
  r = avg(alpha)*avg(h)**2*inner(jump(grad(u),n), jump(grad(v),n))*dS
  
  F = A + r
  
  # Create files for storing results
  file = File("results_%s/u.xdmf" % (fileName))
  
  u = Function(V)
  ffc_options = {"optimize": True, "quadrature_degree": 8}
  problem = LinearVariationalProblem(lhs(F),rhs(F), u, [bc], form_compiler_parameters=ffc_options)
  solver  = LinearVariationalSolver(problem)
  
  
  u.assign(u0)
  u.rename("u", "u")
  
  # Time-stepping
  t = 0.0
  
  file << u
  
  while t < t_end:
  
  	print "t =", t, "end t=", t_end
  
  	# Compute
  	solver.solve()
          plot(u)
  	# Save to file
  	file << u
  
  	# Move to next time step
  	u0.assign(u)
  	t += dt
  

• solved in GLS (Galerkin least squares) formulation

  import time
  import os
  import math
  from dolfin import *
  
  # get file name
  fileName = os.path.splitext(__file__)[0]
  
  parameters['form_compiler']['cpp_optimize'] = True
  parameters['form_compiler']['optimize'] = True
  
  # Parameters
  Pe = Constant(1e10)
  t_end = 10 
  dt = 0.1
  
  # Create mesh and define function space
  mesh = RectangleMesh(0, 0, 1, 1, 40, 40, 'crossed')
  
  # Define function spaces
  V = FunctionSpace(mesh, "CG", 1)
  U = FunctionSpace(mesh, "RT", 1)
  W = MixedFunctionSpace([V,U])
  
  #ic= Expression("((pow(x[0]-0.25,2)+pow(x[1]-0.25,2))<0.2*0.2)?(-25*((pow(x[0]-0.25,2)+pow(x[1]-0.25,2))-0.2*0.2)):(0.0)")
  ic= Expression("((pow(x[0]-0.3,2)+pow(x[1]-0.3,2))<0.2*0.2)?(1.0):(0.0)", domain=mesh)
  
  b = Expression(("-(x[1]-0.5)","(x[0]-0.5)"), domain=mesh)
  
  bc=DirichletBC(W.sub(0),Constant(0.0),DomainBoundary())
  	
  # Define unknown and test function(s)
  w = Function(W)
  
  (u,j) = split(w)
  
  u0 = Function(V)
  
  r1=(1/dt)*(u-u0) + inner(b,j) - (1.0/Pe)*div(j)
  r2= j-grad(u)
  
  LS= inner(r1,r1)*dx + inner(r2,r2)*dx
  
  F = derivative(LS, w)
  J = derivative(F,w)
  
  # Create files for storing results
  file = File("results_%s/u.xdmf" % (fileName))
  
  ffc_options = {"optimize": True, "quadrature_degree": 8}
  problem = NonlinearVariationalProblem(F, w, [bc], J, form_compiler_parameters=ffc_options)
  solver  = NonlinearVariationalSolver(problem)
  
  u0.assign(interpolate(ic,V ))
  
  # Time-stepping
  t = 0.0
  
  (u,j)=w.split()
  u.rename("u", "u")
  file << w
  
  plt=plot(u)
  
  while t < t_end:
  
  	print "t =", t, "end t=", t_end
  
  	# Compute
  	solver.solve()
          plt.plot()
          
          (u,j)=w.split(True)
          u.rename("u", "u")
  	# Save to file
  	file << u
  
  	# Move to next time step
  	u0.assign(u)
  	t += dt