This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

DirichletBC on Stokes Demo

+4 votes

I have the standard 3D Stokes demo running fine. Then, I adapted it to a simpler 2D case, which also worked.

But, when I modify the pressure BC from the entire left edge to only the half part of it, applying no-slip conditions to the rest, some pressure peaks appear near the interface of this two BCs.

It happens due to the degree (2) of velocity FunctionSpace? Which method of DirichletBC should I use? (pointwise, geometric or topological)

Here is my complete code:

from dolfin import *

# Test for PETSc or Epetra
if not has_linear_algebra_backend("PETSc") and not has_linear_algebra_backend("Epetra"):
    info("DOLFIN has not been configured with Trilinos or PETSc. Exiting.")
    exit()

if not has_krylov_solver_preconditioner("amg"):
    info("Sorry, this demo is only available when DOLFIN is compiled with AMG "
     "preconditioner, Hypre or ML.");
    exit()

# Load mesh
mesh = UnitSquareMesh(100, 100) # modified

# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q

# Boundaries
def right(x, on_boundary): return x[0] > (1.0 - DOLFIN_EPS)
def upper_left(x, on_boundary): return x[0] < DOLFIN_EPS and x[1] > 0.5 # modified 
def lower_left(x, on_boundary): return x[0] < DOLFIN_EPS and x[1] < 0.5 # modified 
def top_bottom(x, on_boundary):
    return x[1] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS

# No-slip boundary condition for velocity
noslip = Constant((0.0, 0.0)) # modified
bc0 = DirichletBC(W.sub(0), noslip, top_bottom)
bc3 = DirichletBC(W.sub(0), noslip, lower_left) # modified

# Inflow boundary condition for velocity
inflow = Expression(("-sin(x[1]*pi)", "0.0")) # modified
bc1 = DirichletBC(W.sub(0), inflow, right)

# Boundary condition for pressure at outflow
zero = Constant(0)
bc2 = DirichletBC(W.sub(1), zero, upper_left) # modified

# Collect boundary conditions
bcs = [bc0, bc1, bc2, bc3] # modified

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx

# Form for use in constructing preconditioner matrix
b = inner(grad(u), grad(v))*dx + p*q*dx

# Assemble system
A, bb = assemble_system(a, L, bcs)

# Assemble preconditioner system
P, btmp = assemble_system(b, L, bcs)

# Create Krylov solver and AMG preconditioner
solver = KrylovSolver("tfqmr", "amg")

# Associate operator (A) and preconditioner matrix (P)
solver.set_operators(A, P)

# Solve
U = Function(W)
solver.solve(U.vector(), bb)

# Get sub-functions
u, p = U.split()

# Plot solution
plot(u)
plot(p)
interactive()
asked Apr 17, 2014 by amigoricardo FEniCS Novice (340 points)
edited Apr 29, 2014 by amigoricardo

You should provide complete code to reproduce the problem.

Jan, please take a look on the complete code that I just included.

...