Incompressible Navier-Stokes with GMSH: Found no facets matching domain for BC

mesh
1

Hello everyone,

I am replying the test problem 2 for the incompressible Navier–Stokes equations from the FEniCS tutorial (file ft08_navier_stokes_cylinder.py).
During my short experience with FEniCS, I have solved various PDEs systems by employing GMSH as my Finite-element mesh generator. However, it is the first time I face the issue I am presenting below.

First, to transform GMSH meshes into the correct format, I always effectuate the following steps:

  1. Generate mesh in format .geo followed by .msh one.
  2. Run the following command in Linux terminal to transform .msh format into .xdmf:

meshio-convert Gmshmesh.msh Gmshmesh.xdmf --prune

  1. Run the following command in Linux terminal to eliminate z-coordinate (since I am solving a 2D system):

meshio-convert Gmshmesh.xdmf Gmshmesh.xdmf --prune-z-0

Secondly, this is the Python code for the incompressible Navier–Stokes equations:

from dolfin import Mesh, XDMFFile, MPI, VectorFunctionSpace, FunctionSpace, \
                   DirichletBC, Expression, Constant, TrialFunction, \
                   TestFunction, Function, FacetNormal, sym, nabla_grad, \
                   Identity, dot, inner, div, dx, ds, assemble, lhs, rhs, \
                   TimeSeries, solve, Point, SubDomain, near

from mshr import Rectangle, Circle, generate_mesh

t_final = 5.0             # final time
num_steps = 5000          # number of time steps
dt = t_final / num_steps        # time step size
mu = 0.001                # dynamic viscosity
rho = 1                   # density

"""
# Create mesh by MSHR 
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)
# """

# Reading mesh built by GMSH
mesh = Mesh()
filename = "Gmshmesh.xdmf"
with XDMFFile(MPI.comm_world, filename) as infile:
    infile.read(mesh)

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
walls = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow())
bcu_walls = DirichletBC(V, Constant((0, 0)), walls())
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder())
bcp_outflow = DirichletBC(Q, Constant(0), outflow())
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)

# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
f = Constant((0, 0))
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Time-stepping
t = 0
for n in range(num_steps):
    print('{} out of {}'.format(n, num_steps))

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)

This is the Output error:

*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
0 out of 5000
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
...
...
7 out of 5000
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
Traceback (most recent call last):
  File "Navier_Stokes08.py", line 178, in <module>
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')
  File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 227, in solve
    return dolfin.la.solver.solve(*args)
  File "/usr/lib/python3/dist-packages/dolfin/la/solver.py", line 72, in solve
    return cpp.la.solve(A, x, b, method, preconditioner)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** 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: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

I have reviewed some similar Topics (1, 2, 3, 4) and also rewritten the boundaries definition through Python classes instead of C++ expressions, but I still obtain the same error.

This example runs without any inconvenient when using mshr as the mesh generator (see python code above: lines disabled).

  • Any idea about a possible reason why FEniCS does not recognize certain (GMSH) mesh properties in the Navier–Stokes test problem ??

I apologize for the Post extension; I intend that the error can be replicated or at least well understood. Thank you in advance for your time, and I look forward to hearing from you.

Santiago

2

I am not at a computer right now, so cant tell you why your snippet doesn’t work. However, as you use Gmsh, you are able to mark the Facets for each boundary there. These can be loaded to dolfin with meshio. See Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio for a full code snippet of how to do it with 3D meshes. (Small modifications needed for 2D: i.e. Tetra->triangle, triangle->line, dimension in meshvaluecollection 2->1)