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 (https://fenicsproject.org/pub/tutorial/python/vol1/ft08_navier_stokes_cylinder.py).
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)
set_log_level(1)
# 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))
print(result)
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)
tree_conditions.append(on_tree)
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")
u_n.assign(u_)
p_n.assign(p_)
Full error description:
Traceback (most recent call last):
File "3dstocksminimal.py", line 173, in <module>
solve(A1, u_.vector(), b1, "bicgstab", "default")
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 227, in solve
return dolfin.la.solver.solve(*args)
File "/usr/local/lib/python3.6/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|| = 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,
Konrad
