Treatment of corner nodes with different types of boundary conditions on each side

Hi,

I have a 2D rectangular system with 3 different boundary conditions. The left and right edges have a Dirichlet boundary condition, $ u = g_1 $, and the bottom edge also has a Dirichlet boundary condition, $ u = g_2(x) $. The top boundary is a mixed boundary condition, $ u_x = h_1 $, $ \partial_y u_y = h_2 $.

My understanding is that the corner nodes where the boundaries meet should have the same boundary condition as the left and right edges. Currently, I have the Dirichlet BCs set up using the DirichletBC class and the Neumann part of the mixed BC is in the weak form. However, I'm having problems with the solution near the top corners where the solution is diverging:
screenshot

How does FEniCS treat the nodes at the top corners where the two different types of boundary conditions meet and is there a way to ensure that the corner nodes will have the Dirichlet BC $ u = g_1 $? Is there some special treatment required for cases like this? If so, can you point me towards some literature on this?

Thank you.

ETA: The system is a two-phase flow system where there are two sets of governing equations, one for each phase. I'm using an explicit scheme where I first compute the tentative velocities. I then use the velocities to compute the pressure and then use the pressure to compute the velocities at the current time step. From that, I use the conservation of mass equation to get the phase fractions (alpha): $ \frac{\partial \alpha_i}{\partial t} + \nabla \cdot (\alpha_i u_i) = 0 $. The screenshot shown above is from a plot of one of the alphas. I've included the code in the comment below.

Hi,

Here you go:

from dolfin import *

class u_in_a(Expression):
    def __init__(self, max_vel):
        self.max_vel = max_vel
    def eval(self, values, x):
        values[0] = 0.
        values[1] = self.max_vel*(1-((x[0]*x[0])/(0.1*0.1)))
    def value_shape(self):
        return (2,)

class u_in_b(Expression):
    def __init__(self,max_vel):
        self.max_vel = max_vel
    def eval(self, values, x):
        values[0] = 0.
        values[1] = self.max_vel*(1-((x[0]*x[0])/(0.1*0.1)))
    def value_shape(self):
        return (2,)

class p_rhogh(Expression):
    def __init__(self, rho_init, height):
        self.rho_init = rho_init
        self.height = height
    def eval(self, values, x):
        values[0] = 0. + self.rho_init*9.81*(self.height - x[1])

comm = mpi_comm_world()
mpirank = MPI.rank(comm)

# Load mesh from file
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), "2D.h5", "r")
hdf.read(mesh, "/mesh", False)
boundaries = FacetFunction("size_t", mesh)
hdf.read(boundaries, "/boundaries")

# Define function spaces
# velocities
V1 = VectorFunctionSpace(mesh, "CG", 2)
V2 = VectorFunctionSpace(mesh, "CG", 2)
# alphas
B1 = FunctionSpace(mesh, "CG", 1)
B2 = FunctionSpace(mesh, "CG", 1)
# pressure
Q = FunctionSpace(mesh, "CG", 1)
# mixed space for velocity
V = MixedFunctionSpace([V1, V2])
# mixed space for alphas
B = MixedFunctionSpace([B1, B2])

# define functions
u = TrialFunction(V)
(u_a, u_b) = split(u)
(v_a, v_b) = TestFunctions(V)

alphas = TrialFunction(B)
(alpha_a, alpha_b) = split(alphas)
(b_a, b_b) = TestFunctions(B)

u0 = Function(V)
(u_a0, u_b0) = split(u0)

alphas0 = Function(B)
(alpha_a0, alpha_b0) = split(alphas0)

p = TrialFunction(Q)
q = TestFunction(Q)

p0 = Function(Q)

# time step
timestep = 1.e-5
t_stop = 5.
dt = Constant(timestep)
t = 0.

# define physical properties
rho_a = Constant(1000.)
nu_a = Constant(1e-6)
rho_b = Constant(1.)
nu_b = Constant(1.6e-5)

g = Constant((0., -9.81))

# ICs
zero_vector = Constant((0., 0.))
alpha_a_init = 0.8
alpha_b_init = 0.2

alpha_init = Constant((alpha_a_init, alpha_b_init))
u_init = Constant((0., 0., 0., 0.))

rho_init = alpha_a_init * rho_a + alpha_b_init * rho_b

height = 1.

p_init = p_rhogh(rho_init=rho_init, height=height)

u0.interpolate(u_init)

alphas0.interpolate(alpha_init)

p0.interpolate(p_init)

# create functions for storing the solutions
u1 = Function(V)
(u_a1, u_b1) = split(u1)
u1.interpolate(u_init)

u2 = Function(V)
(u_a2, u_b2) = split(u2)
u2.interpolate(u_init)

alphas1 = Function(B)
(alpha_a1, alpha_b1) = split(alphas1)
alphas1.interpolate(alpha_init)

p1 = Function(Q)
p1.interpolate(p_init)

# Define boundary conditions
# 1 - outlet
# 2 - inlet
# 3 - walls

max_vel_a = 0.02 * alpha_a_init
max_vel_b = 0.02 * alpha_b_init

inflow_a = u_in_a(max_vel_a)
inflow_b = u_in_b(max_vel_b)

noslip1 = DirichletBC(V.sub(0), zero_vector, boundaries, 3)
noslip2 = DirichletBC(V.sub(1), zero_vector, boundaries, 3)

inflow1 = DirichletBC(V.sub(0), inflow_a, boundaries, 2)
inflow2 = DirichletBC(V.sub(1), inflow_b, boundaries, 2)

inflow3 = DirichletBC(B.sub(0), alpha_a_init, boundaries, 2)
inflow4 = DirichletBC(B.sub(1), alpha_b_init, boundaries, 2)

outflow = DirichletBC(Q, 0., boundaries, 1)

outflow11 = DirichletBC(V.sub(0).sub(0), 0., boundaries, 1)
outflow21 = DirichletBC(V.sub(1).sub(0), 0., boundaries, 1)

n = FacetNormal(mesh)
ds = Measure("ds")[boundaries]

I = Identity(u_a.geometric_dimension())

zero_n = Constant((1., 0.))     # this vector will be used in element-wise multiplication to zero out the normal component of the outlet

def epsilon(u_):
    return 0.5*sym(nabla_grad(u_))

def stress(u_, p_, rho_, nu_):
    return 2 * nu_ * epsilon(u_) - nu_ * (2./3.) * I * nabla_div(u_) - p_ * I / rho_

def outflow_stress(u_, p_, rho_, nu_):
    n_stress = dot(n, stress(u_, p_, rho_, nu_))
    # y is normal to the outlet - the normal component of n dot stress = 0
    return elem_mult(n_stress, zero_n)

U_a = 0.5*(u_a + u_a0)
U_b = 0.5*(u_b + u_b0)

# tentative velocity equations
F1a = (1. / dt) * inner(u_a - u_a0, v_a) * dx \
    + inner(dot(u_a0, nabla_grad(u_a0)), v_a) * dx \
    + inner(stress(U_a, p0, rho_a, nu_a), epsilon(v_a)) * dx \
    - inner(outflow_stress(U_a, p0, rho_a, nu_a), v_a) * ds(1) \
    - inner(g, v_a) * dx \
    - (1. / rho_a) * inner(0.05 * (u_b - u_a) * alpha_b0 / alpha_a0, v_a) * dx

F1b = (1. / dt) * inner(u_b - u_b0, v_b) * dx \
    + inner(dot(u_b0, nabla_grad(u_b0)), v_b) * dx \
    + inner(stress(U_b, p0, rho_b, nu_b), epsilon(v_b)) * dx \
    - inner(outflow_stress(U_b, p0, rho_b, nu_b), v_b) * ds(1) \
    - inner(g, v_b) * dx \
    + (1. / rho_b) * inner(0.05 * (u_b - u_a), v_b) * dx

F1 = F1a + F1b

a1 = lhs(F1)
L1 = rhs(F1)

# pressure update equation
F2 = (rho_a / dt) * nabla_div(alpha_a0 * u_a1) * q * dx \
    + (rho_b / dt) * nabla_div(alpha_b0 * u_b1) * q * dx \
    + inner(nabla_grad(p - p0), nabla_grad(q)) * dx

a2 = lhs(F2)
L2 = rhs(F2)

# velocity update
F3a = (1. / dt) * inner(u_a - u_a1, v_a) * dx \
    + (1. / rho_a) * inner(nabla_grad(p1), v_a) * dx

F3b = (1. / dt) * inner(u_b - u_b1, v_b) * dx \
    + (1. / rho_b) * inner(nabla_grad(p1), v_b) * dx

F3 = F3a + F3b

a3 = lhs(F3)
L3 = rhs(F3)

# alpha update
F4a = (1. / dt) * (alpha_a - alpha_a0) * b_a * dx \
    + nabla_div(alpha_a * u_a2) * b_a * dx

F4b = (1. - alpha_a - alpha_b) * b_b * dx

F4 = F4a + F4b

a4 = lhs(F4)
L4 = rhs(F4)

problems = [
    LinearVariationalProblem(a1, L1, u1, [noslip1, noslip2, inflow1, inflow2, outflow11, outflow21]),
    LinearVariationalProblem(a2, L2, p1, [outflow]),
    LinearVariationalProblem(a3, L3, u2, [noslip1, noslip2, inflow1, inflow2, outflow11, outflow21]),
    LinearVariationalProblem(a4, L4, alphas1, [inflow3, inflow4]),
]
solvers = [
    LinearVariationalSolver(problems[0]),
    LinearVariationalSolver(problems[1]),
    LinearVariationalSolver(problems[2]),
    LinearVariationalSolver(problems[3]),
]


def solve_step():
    begin("Computing tentative velocity")
    solvers[0].solve()
    end()

    # Pressure correction
    begin("Computing pressure update")
    solvers[1].solve()
    end()

    # Velocity correction
    begin("Computing velocity update")
    solvers[2].solve()
    end()

    # alpha update
    begin("Computing alpha update")
    solvers[3].solve()
    end()

    u0.assign(u2)
    p0.assign(p1)
    alphas0.assign(alphas1)


while t < t_stop + DOLFIN_EPS:

    if mpirank == 0:
        print "t = ", t

    if mpirank == 0:
        print common.memory_usage()

    t += timestep

    solve_step()

Thank you.