building on the instationary advection-diffusion demo I used the Galerkin approach with SUPG stabilisation terms in order to solve the same scalar transport equation for stationary conditions. The solver works well as long as the diffusive term is small (i.e. |c/velocity| <~ 0.05). When all boundaries are periodic (which I need) the solution looks even worse. Is there a way to improve the solution?

Here is my MWE (runs on Linux for the latest stable version of Fenics to the date)

```
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#Dimensions
Lx = 1
Ly = 1
Nx = 100
class PeriodicBoundaries(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the two corners (0, 1) and (1, 0)
return bool((near(x[0], 0) or near(x[1], 0)) and
(not ((near(x[0], 0) and near(x[1], Ly)) or
(near(x[0], Lx) and near(x[1], 0)))) and on_boundary)
def map(self, x, y):
if near(x[0], Lx) and near(x[1], Ly):
y[0] = x[0] - Lx
y[1] = x[1] - Ly
elif near(x[0], Lx):
y[0] = x[0] - Lx
y[1] = x[1]
elif near(x[1], Ly):
y[0] = x[0]
y[1] = x[1] - Ly
else:
y[0] = -1000
y[1] = -1000
# Sub domain for Periodic boundary condition
class PeriodicBoundary_y(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return bool(x[1] < DOLFIN_EPS and x[1] > -DOLFIN_EPS and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
y[0] = x[0]
y[1] = x[1] - Ly
# Define mesh
mesh = UnitSquareMesh(Nx, Nx)
# Build function space
#~ Q = FunctionSpace(mesh, "Lagrange", 1, constrained_domain = PeriodicBoundaries())
#~ Q = FunctionSpace(mesh, "Lagrange", 1, constrained_domain = PeriodicBoundary_y())
Q = FunctionSpace(mesh, "Lagrange", 1)
# Define velocity as vector
velocity = as_vector([1.0, 0.0])
#~ velocity = as_vector([0.0, 1.0])
c = Constant(0.05)
# Test and trial functions
u, v = TrialFunction(Q), TestFunction(Q)
# Manufactured right hand side
x = SpatialCoordinate(mesh)
R = 0.303825*Lx # Radius of the source in the center of the mesh
r1 = np.power( np.power( x[0]- 0.5 , 2 ) + np.power( x[1]- 0.5 , 2 ), 0.5 )
f = conditional(gt(r1, R), 0, 1.0) # circular source
# Residual
r = dot(velocity, grad(u)) - c*div(grad(u)) - f
# Galerkin variational problem
F = v*dot(velocity, grad(u))*dx + c*dot(grad(v), grad(u))*dx - f*v*dx
# Add SUPG stabilisation terms
vnorm = sqrt(dot(velocity, velocity))
h = CellDiameter(mesh)
delta = h/(2.0*vnorm)
F += delta*dot(velocity, grad(v))*r*dx
# Create bilinear and linear forms
a = lhs(F)
L = rhs(F)
A, b = assemble_system(a, L)
u = Function(Q)
bc = []
# Solve
u = Function(Q)
solve(a == L, u, bc)
plot(u)
plt.show()
```