Failed to converge in incompressible Navier-Stokes of example from lecture 9, 2013-11-logg

Question

I am trying to do the challenge from lecture 9 of the course in Nov-2013 by Logg (http://fenicsproject.org/pub/course/lectures/). The challenge is to implement the incompressible Navier-Stokes on the dolphin mesh, with 1kPa at inflow, and 0kPa at outflow. On the rest of the boundaries, the velocity is set to 0. My approach was to simply start from the sample code in the tutorial for the Lshape tube (tutorial here), making only the few following changes:

# Load mesh from file
mesh = Mesh("../dolfin_fine.xml.gz")
sub_domains = MeshFunction("size_t", mesh, "../dolfin_fine_subdomains.xml.gz")

# Set parameter values
dt = 0.0005
T = 0.1
nu = 0.001002/1000

# Define boundary conditions
noslip = DirichletBC(V, Constant((0,0)), sub_domains, 0)

inflow = DirichletBC(Q, Constant(1000), sub_domains, 1)

outflow = DirichletBC(Q, Constant(0), sub_domains, 2)

The mesh and its boundary definitions can be downloaded here.
The solution fails to converge, with the following message:

*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = inf).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: unknown

Could someone tell me what's wrong with my implementation? Thank you in advance.

Answer 1

I just found the mistake: the code in the tutorial was actually implementing the Navier-Stokes equation with unity density, therefore it didn't work for the problem I was trying to solve (water properties). After the small correction, the code worked.