Boundary condition not enforced

1

Hi,
I am solving for the pressure of a non-Newtonian fluid model and the boundary condition is not enforced. We specify the pressure on one side of the domain and solve, but the recovered solution is zero pretty much everywhere. I’ve tried iterative solvers with nonzero initial guesses but don’t see a difference. Maybe I’m missing something really obvious?

Thanks to anyone who gives this a look. I’ve been fiddling with this for a long time and am up for tying anything!

from dolfin import *

mesh = RectangleMesh(Point(-4, 0.0), Point(4, 1.0), 32, 32, "crossed")

# refine the boundary
for j in range(6):
	mesh.init(1, 2)
	markers = CellFunction("bool", mesh, False)

	for c in cells(mesh):
		for f in facets(c):
			if f.exterior():
				markers[c] = True
				break

	mesh = refine(mesh, markers)

# create function spaces
V = VectorFunctionSpace(mesh, "Lagrange", 4)
P = FunctionSpace(mesh, "Lagrange", 3)

boundary_exp_u = Expression(("6*x[1]*(1 - x[1])", "0.0"))
boundary_exp_p = Expression("-9.6*x[0]")

def p_boundary(x, on_boundary):
	return on_boundary and x[0] > 4 - 3*DOLFIN_EPS

bc_u = DirichletBC(V, boundary_exp_u, "on_boundary")
bc_p = DirichletBC(P, boundary_exp_p, p_boundary)

# exact solutions
u = project(boundary_exp_u, V)
pi = project(boundary_exp_p, P)

# pressure solve
p = TrialFunction(P)
q = TestFunction(P)

a = (p*(q + 5.0*dot(u, grad(q))))*dx
L = pi*q*dx

p = Function(P)	
problem = LinearVariationalProblem(a, L, p, bc_p)
solver = LinearVariationalSolver(problem)
solver.solve()

plot(p, interactive=True)
2

When I open the solution in Paraview, it looks like the BC is enforced correctly (although I’m using 2018.1, and it looks like you’re using an older version). The rest of the solution is really bad, though, because the variational problem is unstable. It is reaction–advection in conservative form, but missing a boundary term. The following change to the bilinear form leads to a stable formulation, and a solution where it is more visually clear that the BC is being enforced:

# Inflow stabilization
n = FacetNormal(mesh)
a += dot(-5.0*u,n)*p*q*ds # (assuming Dirichlet BCs at inflow, div(u)=0)

# General case, if you have inflow Neumann boundaries:
#from ufl import Max # (Not necessary in older versions)
#a += Max(dot(-5.0*u,n),Constant(0.0))*p*q*ds
3

Thanks for your advice, the solution looks much better with the stabilization. I’m not sure exactly what was going on with the boundary but I was using an iterative solver, which seems to be no good for this problem. Thanks!