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