Improve convergence/stability in stationary scalar transport equation

variational formulation
1

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()
2

The problem becomes ill-posed in two ways if the domain is fully periodic. First, the fully-periodic problem only makes sense if the integral of f is zero (so u is being produced and consumed at the same rate, allowing some finite steady-state solution). Second, the solution u is only defined up to a constant. See the attached modification of the MWE, in which the integral of u is constrained to equal some value, and f is corrected to have zero integral:

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
VE = FiniteElement("Lagrange",mesh.ufl_cell(),1)
# Element for Lagrange multiplier to enforce integral constraint:
LE = FiniteElement("Real",mesh.ufl_cell(),0)
Q = FunctionSpace(mesh,MixedElement([VE,LE]),
                  constrained_domain = PeriodicBoundaries())
#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)
u = U[0]
lam = U[1]
v = V[0]
rho = V[1]

# 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

# Fix `f` so the problem is well-posed:
ftotal = assemble(f*dx(domain=mesh))
f = f - ftotal

# Residual
r = dot(velocity, grad(u)) - c*div(grad(u)) - f

# Choose some value that the total integral of `u` should be equal to:
utotal = Constant(1.0)

# Galerkin variational problem
F = v*dot(velocity, grad(u))*dx + c*dot(grad(v), grad(u))*dx - f*v*dx \
    + lam*v*dx - rho*(u-utotal)*dx # Enforce value for integral of `u`

# 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)
u,lam = U.split()
plot(u)
plt.show()
3

Thanks for the support. It works.