Convergence error for 3D Navier-Stokes after changing mesh

Hi all,

I’m using FEniCS to simulate the 3D Navier Stokes equation. The geometry are plotted in the picture below. The code I use is a modified example from tutorial (
The main goal of my project is carry out a forest fire simulation, and thoroughly examine how the combustion products propagate, especially PM dust (I am working on a sensor to detect wild fires).

In the first stage, I want to find pressure and velocity by solving 3D NS equation, in a “forest” domain (3DBox minus simple figures resembling trees).

Problem: The solver diverges after few timesteps. I have tried increasing mesh resolution and shortening time step with no luck. The same program on a simple mesh (only 3D BOX without cutting cylinders and spheres from domain) works perfectly fine.

Things I checked or tried to change:

  • I checked the coverage of the boundary conditions, no warning Found no facets matching domain for boundary condition

  • I increased segments parameter in the constructor of Cylinder and Sphere class

  • I add to generate_mesh function tetgen parameter

  • I have tried increasing mesh resolution and shortening time step. The only difference was the appearance of a error after different number of time steps

  • I suspected that there may be a contradiction between boundary conditions at the border of inflow and walls. I changed inflow profile function that determines inflow velocity to be 0 on edges.
    (1., 0., 0.) -> (‘fabs(x[1])*fabs(x[1]-4.0)*fabs(x[2])*fabs(x[2]-5.0)’, ‘0’, ‘0’). But if this was the cause of the error, the version with a simple domain (only 3D Box) should also appear error. That’s why I suspect it’s a mesh issue

  • I tried to run the program single-threaded and in parallel with mpirun

  • I also try with direct solver, I remove preconditioner setings and solve function additional parameters.

Unfortunately any of these changes do not eliminated the problem.

Below is a picture of mesh created in jupyter notebook by IPython.display.HTML(X3DOM.html(mesh))


Code of my simulation:

from fenics import *
from mshr import *

T = 1  # final time
num_steps = int(200)  # number of time steps
dt = T / num_steps  # time step size
mu = 0.001  # dynamic viscosity
rho = 1  # density

def generate_tree(x, y, r=0.15, R=0.3, h=1.5):
    cylinder_h = h - R
    cylinder = Cylinder(Point(x, y, cylinder_h), Point(x, y, 0.0), r, r)
    treetop = Sphere(Point(x, y, cylinder_h), R)
    tree = cylinder + treetop
    return tree

box = Box(Point(0, 0, 0), Point(7, 4, 5))

tree1_1 = generate_tree(1.5, 0.75, r=0.15, R=0.5, h=3)
tree1_2 = generate_tree(1.5, 2, r=0.15, R=0.5, h=3)
tree1_3 = generate_tree(1.5, 3.25, r=0.15, R=0.5, h=3)

tree2_1 = generate_tree(2.75, 0.75, r=0.15, R=0.5, h=3)
tree2_2 = generate_tree(2.75, 2, r=0.15, R=0.5, h=3)
tree2_3 = generate_tree(2.75, 3.25, r=0.15, R=0.5, h=3)

tree3_1 = generate_tree(4., 0.75, r=0.15, R=0.5, h=3)
tree3_2 = generate_tree(4., 2, r=0.15, R=0.5, h=3)
tree3_3 = generate_tree(4., 3.25, r=0.15, R=0.5, h=3)

tree4_1 = generate_tree(5.25, 0.75, r=0.15, R=0.5, h=3)
tree4_2 = generate_tree(5.25, 2, r=0.15, R=0.5, h=3)
tree4_3 = generate_tree(5.25, 3.25, r=0.15, R=0.5, h=3)

forest = tree1_1 + tree1_2 + tree1_3 + tree2_1 + tree2_2 + tree2_3 + tree3_1 + tree3_2 + tree3_3 + tree4_1 + tree4_2 + tree4_3
domain = box - forest
mesh = generate_mesh(domain, 40)

# Define function spaces for pressure and velocity
V = VectorFunctionSpace(mesh, "Lagrange", degree=2, dim=3)
Q = FunctionSpace(mesh, "Lagrange", 1)

class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (abs(x[0] - 0.) < DOLFIN_EPS)

class Outflow(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (abs(x[0] - 7.) < DOLFIN_EPS)

class Walls(SubDomain):
    def inside(self, x, on_boundary):
        result = on_boundary and ((abs(x[1] - 0.) < DOLFIN_EPS) or (abs(x[1] - 4.) < DOLFIN_EPS) or
                                  (abs(x[2] - 0.) < DOLFIN_EPS) or (abs(x[2] - 5.) < DOLFIN_EPS))
        return result

class Forest(SubDomain):
    def inside(self, x, on_boundary):
        r = 0.15
        R = 0.5
        h = 3

        x_l = [1.5, 2.75, 4., 5.25]
        y_l = [0.75, 2., 3.25]
        tree_conditions = []
        for xt in x_l:
            for yt in y_l:
                on_treetop = (abs(x[0] - xt) < DOLFIN_EPS + R) and (abs(x[1] - yt) < DOLFIN_EPS + R) and (
                        abs(x[2] - (h - R)) < DOLFIN_EPS + R)
                on_treetrunk = (abs(x[0] - xt) < DOLFIN_EPS + r) and (abs(x[1] - yt) < DOLFIN_EPS + r) and (
                        x[2] < DOLFIN_EPS + (h - R))
                on_tree = on_boundary and (on_treetop or on_treetrunk)
        result = any(tree_conditions)
        return result

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)

Inflow().mark(boundaries, 1)
Outflow().mark(boundaries, 2)
Walls().mark(boundaries, 3)
Forest().mark(boundaries, 4)

inflow_profile = ('1.', '0', '0')

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=3), boundaries, 1)
bcp_outflow = DirichletBC(Q, Constant(0.), boundaries, 2)
bcu_walls = DirichletBC(V, Constant((0., 0., 0.)), boundaries, 3)
bcu_forest = DirichletBC(V, Constant((0, 0, 0)), boundaries, 4)
bcu = [bcu_inflow, bcu_walls, bcu_forest]
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, 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)

# Use amg preconditioner if available
prec = "petsc_amg" if has_krylov_solver_preconditioner("amg") else "default"

# Use nonzero guesses - essential for CG with non-symmetric BC
parameters['krylov_solver']['nonzero_initial_guess'] = True

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

t = 0
for n in range(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", "default")

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

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, "bicgstab", "default")


Full error description:

Traceback (most recent call last):
  File "", line 173, in <module>
    solve(A1, u_.vector(), b1, "bicgstab", "default")
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/", line 227, in solve
  File "/usr/local/lib/python3.6/dist-packages/dolfin/la/", line 72, in solve
    return, x, b, method, preconditioner)

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** 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|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

I use docker on ubuntu 19.10.

I would be very grateful if you could help me solve this problem, I have tried to solve this for the last few days, unfortunately without success. Due to this error, the further stage of my research is entirely blocked

Best Regards,

Hello all.
Can somebody help me with this problem? It really blocked my entire project.

Best Regards,

First step would be to test the flow solver with only one tree, instead of a forest. You could Also try to simplify the tree (make it a rectangular block) to ensure that you are setting the correct boundary conditions on the obstacle.

Thank you very much,
I’ll try this way

I investigate these problem more deeply according to your suggestion, and the results are quite surprising to me.

For the most simplifying example, 1 tree as a rectangular block:

forest = Box(Point(2, 2, 0), Point(3, 3, 2.5))

and the following boundary conditions:

class Forest(SubDomain):
    def inside(self, x, on_boundary):
        result = on_boundary\
        and x[0] > 2.0 - DOLFIN_EPS and x[0] < 3.0 + DOLFIN_EPS\
        and x[1] > 2.0 - DOLFIN_EPS and x[1] < 3.0 + DOLFIN_EPS\
        and x[2] > 0.0 - DOLFIN_EPS and x[2] < 2.5 + DOLFIN_EPS
        return result

solver works fine.

Also if I use only one component of my tree model (cylinder or sphere), everything works correctly (with the same boundary conditions):

def generate_tree(x, y, r=0.15, R=0.3, h=1.5):
    cylinder_h = h - R
    cylinder = Cylinder(Point(x, y, cylinder_h), Point(x, y, 0.0), r, r)
    return cylinder

def generate_tree(x, y, r=0.15, R=0.3, h=1.5):
    cylinder_h = h - R
    treetop = Sphere(Point(x, y, cylinder_h), R)
    return treetop

forest = generate_tree(2.5, 2.5, r=0.2, R=0.5, h=2.5)

Only after applying cylinder and sphere union, solver error appear.

cylinder = Cylinder(Point(x, y, cylinder_h), Point(x, y, 0.0), r, r)
treetop = Sphere(Point(x, y, cylinder_h), R)
tree = cylinder + treetop
return tree

Do you have any ideas on how this behavior might be caused?

Best Regards,

Have you tried to visualize the mesh function used to mark the boundaries in for instance paraview? It might be that you are not marking the boundaries correctly. In general, for 3D geometries, i tend to use Gmsh with the opencascade backend to Get the correct boundary markers.

I saved the mesh function with marked boundaries to .pvd format and opened it in paraview. After filtering by value (to get only surface of the tree) I received a coherent grid, without gaps. After adding walls boundaries I get 1 surface, without any empty spaces.

Best Regards,

The mesh quality in the area close to the intersection between the ball and cylinder seems quite bad in my opinion.

1 Like

Seconding @dokken the mesh is extremely poor quality in the eyeball norm. mshr is no longer supported (and I’ve never liked CGAL for generating high quality FEM meshes either). Look into using gmsh or commerical software for such intricate geometries. Generation of quality meshes is not a trivial task.

1 Like

Hi [Konrad_Zuchniak],
Could you share more information on the physics of the system that you are trying to simulate?
This sounds like highly complex problem, assuming you want to also do combustion (multi-phase reactive flow). How are you intending to do this?