Thank you for spotting my mistakes! I have implemented the suggested corrections and added the no-slip boundary conditions. Now I get the following exception: [scfw@comp-pc3237 CH_advection]$ python navier_stokes_CH.py 0.0005 pyop2:INFO Solving linear variational problem... Traceback (most recent call last): File "navier_stokes_CH.py", line 103, in <module> solve(a == L, u, bcs=[bc0, bc1, bc2, bc3]) # 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 147, 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: DIVERGED_NANORINF ---- New code: from __future__ import division from firedrake import * # Model parameters dt = 5.0e-04 # time step Cn=0.01 # Cahn We=0.45 # Weber Pe=4.5 # Peclet Re=100 # Reynolds gravity = False # switch for gravity Fr2=0.1 # Froude number squared beta = 2 # stabilising term rho1 = Constant(1) rho2 = Constant(1) drho_hat = Constant(0) drho_hat.assign(0.5*(1 - (rho2/rho1))) # d rho_hat / d phi alpha = (rho2-rho1)/(rho2+rho1) # Create mesh and define function spaces mesh = UnitSquareMesh(200, 200) P1 = FunctionSpace(mesh, "Lagrange", 1) P2 = VectorFunctionSpace(mesh, "Lagrange", 2) # Taylor-Hood elements for velocity and pressure, P1-P1 elements for phi and mu W = P2*P1*P1*P1 # Dirichlet boundary condition for velocity noslip = Constant((0, 0)) bc0 = DirichletBC(W.sub(0), noslip, 1) bc1 = DirichletBC(W.sub(0), noslip, 2) bc2 = DirichletBC(W.sub(0), noslip, 3) bc3 = DirichletBC(W.sub(0), noslip, 4) # Define test functions chi, theta, omega, psi = TestFunctions(W) # Define functions u = Function(W) # current solution u0 = Function(W) # solution from previous converged step # Split mixed functions v, p, phi, mu = split(u) v0, p0, phi0, mu0 = split(u0) # Compute the chemical potential df/dphi phi = variable(phi) f = 0.25*(phi**2 - 1)**2 dfdphi = diff(f, phi) # create a functions for the two definitions of density rho_hat = Function(P1) rho = Function(P1) ex, ey = unit_vectors(2) # Weak statement of the equations a1 = (1/dt)*omega*phi*dx + omega*div(v*phi0)*dx + (1/Pe)*inner(grad(omega),grad(mu))*dx a2 = psi*mu*dx - (beta/Cn)*phi*psi*dx - Cn*inner(grad(phi),grad(psi))*dx - (1/rho_hat)*drho_hat*psi*p*dx a3 = rho*(1/dt)*inner(chi, v)*dx + rho*inner(dot(v0,nabla_grad(v)), chi)*dx +(1/We)*(-p*div(chi) + inner(chi, phi*grad(mu))-(1/rho_hat)*drho_hat*p*inner(chi, grad(phi0)) )*dx + (1/Re)*inner((nabla_grad(v)[i,j] + nabla_grad(v)[j,i]), nabla_grad(chi)[i,j] )*dx - (2/(3*Re))*div(v)*div(chi)*dx a4 = theta*div(v)*dx + (alpha/Pe)*inner(grad(theta), grad(mu))*dx a = derivative(a1 + a2 + a3 + a4, u) L1 = (1/dt)*omega*phi0*dx L2 = - (beta/Cn)*phi0*psi*dx + (1/Cn)*dfdphi*psi*dx L3 = rho*(1/dt)*inner(chi, v0)*dx if gravity: L3 += -(1/Fr2)*rho_hat*inner(chi, ey)*dx L = L1 + L2 + L3 v0, p0, phi0, mu0 = u0.split() # Initial conditions for phi phi_ini = Expression('''1 - tanh( ( sqrt( pow(x[0]-x0, 2)+pow(x[1]-y0, 2) ) - r0 )/C*sqrt(2) )- tanh( ( sqrt( pow(x[0]-x1, 2)+pow(x[1]-y1, 2) ) - r1 )/C*sqrt(2) )''', x0=0.4,y0=0.5,r0 = 0.25,x1=0.78,y1=0.5,r1 = 0.1,C=0.01 ) phi0.interpolate(phi_ini) # Initial conditions for velocity vel_ini = Expression(("0.0", "0.0")) v0.interpolate(vel_ini) # initialise density rho.interpolate(0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0))) outfile = File("./Results_NSCH/nsch_res.pvd") t = 0 tend = 1.0 v, p, phi, mu = u0.split() v.rename("Velocity") p.rename("Pressure") phi.rename("Order parameter") mu.rename("Chemical potential") outfile.write(v, p, phi, mu, time=t) while t <= tend: t+=dt ; print t # compute density using algebraic expression rho_hat = 0.5*( (1+phi) + (rho2/rho1)*(1-phi)) rho_hat = 0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0)) # Solve the system of equations solve(a == L, u, bcs=[bc0, bc1, bc2, bc3]) # # update rho rho -= dt*div(rho*v0) u0.assign(u) outfile.write(v, p, phi, mu, time=t) ----