Hello again, 

I'm having some problems trying to implement a symmetry axis boundary condition for a 2d flow (plane poiseuille). My weak form and approach is described in the pdf.

When I try to impose no normal flow through the symmetry axis by setting

bc2 = DirichletBC(W.sub(0).sub(0), 0, 1)

I get the following exception

***

Traceback (most recent call last):

  File "ns_plane_poiseuille.py", line 104, in <module>

    solve(F_u == 0, u, bcs = [bc0,bc1,bc2],solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'} ) # 

  File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve

    _solve_varproblem(*args, **kwargs)

  File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 164, in _solve_varproblem

    solver.solve()

  File "<decorator-gen-295>", line 2, in solve

  File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper

    return f(*args, **kwargs)

  File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", line 190, in solve

    solving_utils.check_snes_convergence(self.snes)

  File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", line 62, in check_snes_convergence

    %s""" % (snes.getIterationNumber(), msg))

RuntimeError: Nonlinear solve failed to converge after 0 nonlinear iterations.

Reason:

   Inner linear solve failed to converge after 0 iterations with reason: unknown reason (petsc4py enum incomplete?)

***

If I run the code without that boundary condition, I get unwanted flow through that symmetry axis....

Any advice on setting the symmetry bc would be greatly appreciated

Full code:

from __future__ import division

from firedrake import *

import numpy as np

class Configuration(object):

    def __init__(self, **kwargs):

        for name, value in kwargs.iteritems():

            self.__setattr__(name, value)

    def __setattr__(self, name, value):

        """Cause setting an unknown attribute to be an error"""

        if not hasattr(self, name):

            raise AttributeError("'%s' object has no attribute '%s'" % (type(self).__name__, name))

        object.__setattr__(self, name, value)

class MyParameters(Configuration):

# Model parameters

dt = 1.0e-02 # time step

Cn=0.03 # Cahn

We=1 # Weber

Pe=20000 # Peclet

Re=1 # Reynolds

gravity = False # switch for gravity 

Fr2=0.1 # Froude number squared

beta = 2 # stabilising term

rho1 = Constant(1.0)

rho2 = Constant(1.0)

drho_hat = 0.5*(1 - (rho2/rho1))

alpha = (rho2-rho1)/(rho2+rho1)

class phiBC(DirichletBC):

def __init__(self, V, value, sid):

super(phiBC,self).__init__(V,value,sid)

self.nodes = np.arange(V.dof_dset.total_size, dtype=np.int32)

myparameters = MyParameters()

# Create mesh and define function spaces

m=64

mesh = UnitSquareMesh(m, m)

P1 = FunctionSpace(mesh, "Lagrange", 1)

P2 = VectorFunctionSpace(mesh, "Lagrange", 2)

n = FacetNormal(mesh)

# Create mixed space

W = P2*P1

# Define functions

u   = Function(W) # 

u0  = Function(W) # lagged solution

# Dirichlet boundary condition for velocity

noslip = Constant((0, 0))

inlet = Function(P2).interpolate(Expression(("0.0","0.25*(x[0] + 1)*(x[0] - 1)")))

bc0 = DirichletBC(W.sub(0), inlet, 4)

bc1 = DirichletBC(W.sub(0), noslip, 2)

bc2 = DirichletBC(W.sub(0).sub(0), 0, 1)

# Initial conditions

v0, p0 = u0.split()

vel_ini = Expression(("0.0", "0.0")) # Initial conditions for velocity

v0.interpolate(vel_ini)

def F(u, u0, rho, rho0, rho_hat, W,n, myparameters):

# access parameters

Re=myparameters.Re

# Define test functions

chi, theta = TestFunctions(W)

v, p  = split(u)

v0, p0 = split(u0)

f3 = inner(dot(v0,nabla_grad(v)), chi)*dx \

+ (-p*div(chi) )*dx \

+ (1/Re)*inner(nabla_grad(v), nabla_grad(chi) )*dx \

+ (1/Re)*inner(nabla_grad(v)[j,i], nabla_grad(chi)[i,j])*dx - (2/(3*Re))*div(v)*div(chi)*dx \

- (1/Re)*inner(inner(chi[i], nabla_grad(v)[i,j]), n[j])*(ds(1)+ds(4)+ds(3)) \

f4 = theta*div(v)*dx 

F = f3 + f4

return F

v, p = u0.split()

v.rename("Velocity")

p.rename("Pressure")

outfile = File("./Results/nsch_res.pvd")

parameters["matnest"] = False

exact = Function(P2).interpolate(inlet)

rho1 = myparameters.rho1

rho2 = myparameters.rho2

F_u = F(u,u, rho1, rho1, rho1, W, n, myparameters)

t = 0.0

# Solve the system of equations

solve(F_u == 0, u, bcs = [bc0,bc1,bc2],solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'} ) # 

u0.assign(u)

v, p= u0.split()

outfile.write(v, p, time=t)

print "L2norm = ", sqrt(assemble(inner(v-exact,v-exact)*dx))