Gmsh generated mesh - solver doesnt converge -code included

Question

Hi, I am solving a simple 2d Navier Stokes problem. I import a mesh exported from gmsh as shown,

My solver crashes with error: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = -nan).

The code,

 from dolfin import *
import numpy as np
set_log_level(WARNING)
nx=1250.
ny=750.
class NoslipBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary #and (near(x[1],0) or near(x[1],ny))

class InflowBoundary(SubDomain):
    def inside(self, x, on_boundary):        
        return near(x[0],0) and on_boundary
class OutflowBoundary(SubDomain):
    def inside(self, x, on_boundary):        
        return near(x[0],nx) and on_boundary


def flowin():  

  mesh = Mesh("square.xml")
  print Mesh.num_cells(mesh)
  h = mesh.hmin()
  f = Constant((0., 0.))         # force
  nu = Constant(1./10)           # kinematic viscosity  
  dt = 0.2*h/0.3               # time step CFL with 1 = max. velocity
  k = Constant(dt)               # time step 
  T = 200.0                       # total simulation time
  u0 = Constant((0., 0.))        # initial velocity
  p0 = Constant(0.)             # initial pressure


  V = VectorFunctionSpace(mesh, 'CG', 2)
  Q = FunctionSpace(mesh, 'CG', 1)

  u = TrialFunction(V)
  v = TestFunction(V)
  p = TrialFunction(Q)
  q = TestFunction(Q)

  u0 = interpolate(u0, V)  
  p0 = interpolate(p0, Q)
  us = Function(V)
  u1 = Function(V)
  p1 = Function(Q)

  # tentative velocity
  F0 = (1./k)*inner(u - u0, v)*dx + inner(dot(grad(u0), u0), v)*dx\
       + nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
  a0, L0 = system(F0)

  # projection
  F1 = inner(grad(p), grad(q))*dx + (1./k)*q*div(us)*dx
  a1, L1 = system(F1)

  # finalize
  F2 = (1./k)*inner(u - us, v)*dx + inner(grad(p1), v)*dx 
  a2, L2 = system(F2)

  # bcs
  b_v = DirichletBC(V, Constant((0., 0.0)), NoslipBoundary())  
  b_v1 = DirichletBC(V, Constant((0.01,0)), InflowBoundary())

  b_p1 = DirichletBC(Q, Constant(0.),  OutflowBoundary())

  bcs_v = [b_v,b_v1]
  bcs_p = [b_p1]

  A0 = assemble(a0)
  A1 = assemble(a1)
  A2 = assemble(a2)
  prec = "amg" if has_krylov_solver_preconditioner("amg") else "default" 
  solver02 = KrylovSolver('gmres', 'ilu')
  solver1 = KrylovSolver('gmres', 'petsc_amg')
  solver3 = KrylovSolver('gmres', 'ilu')
  t = 0
  count=0
  while t < T:

    count=count+1
    b = assemble(L0)
    [bc.apply(A0, b) for bc in bcs_v]
    solver02.solve(A0, us.vector(), b)

    b = assemble(L1)
    [bc.apply(A1, b) for bc in bcs_p]
    solver1.solve(A1, p1.vector(), b)


    b = assemble(L2)
    [bc.apply(A2, b) for bc in bcs_v]
    solver02.solve(A2, u1.vector(), b)
    u0.assign(u1)

    t += dt
    plot(u0, title="Velocity", rescale=True, interactive=True)     
  return u0, p1

#----------------------------------------------------------------------

if __name__ == '__main__':
  u, p = flowin()

The code works fine if I just have a rectangular domain without the hole at the center. I have tried reducing the inlet velocity and time step, but no change. Could someone help me find the mistake with the code or solver? Thanks in advance.

Answer 1

Your boundary conditions are wrong, you are enforcing noslip on the outlet. Otherwise your code works for me on a similar mesh generated with mshr.

Comments

Thank you, I will check my outlet boundary conditions.

---

You could try

class NoslipBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and not (near(x[0],0) or near(x[0],nx))

---

Now, it works perfectly. Thank you