Apologies, here's the bitbucket repo https://bitbucket.org/fryderyk216/navierstokes the relevant file is called ns_plane_poiseuille . py Thanks, Fryderyk ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of David Ham <David.Ham@imperial.ac.uk> Sent: 25 April 2016 11:11:46 To: firedrake Subject: Re: [firedrake] Symmetry axis boundary condition Hi Fryderyk, Pasting Python code into emails is always unsafe. In this case, the indentation has been mucked up, which changes the meaning of the code. Can you please share a git repo (eg on bitbucket or github) with the code so I can try running it. Regards, David On Sun, 24 Apr 2016 at 13:47 Fryderyk Wilczynski <scfw@leeds.ac.uk<mailto:scfw@leeds.ac.uk>> wrote: 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))