# 2D buoyancy-driven flow (Boussinesq) through inlet patch at the top wall with circular obstacle in the room

**URL:** <https://fenicsproject.discourse.group/t/2d-buoyancy-driven-flow-boussinesq-through-inlet-patch-at-the-top-wall-with-circular-obstacle-in-the-room/18210>\
**Category:** Errors\
**Created:** [September 15, 2025, 6:34am UTC](https://fenicsproject.discourse.group/t/2d-buoyancy-driven-flow-boussinesq-through-inlet-patch-at-the-top-wall-with-circular-obstacle-in-the-room/18210 "2025-09-15T06:34:23Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Abrar\_Nhavakar](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/abrar_nhavakar/32/8064_2.png) [@Abrar\_Nhavakar](https://fenicsproject.discourse.group/u/Abrar_Nhavakar)\
**Post date:** [September 15, 2025, 6:34am UTC](https://fenicsproject.discourse.group/t/2d-buoyancy-driven-flow-boussinesq-through-inlet-patch-at-the-top-wall-with-circular-obstacle-in-the-room/18210/1 "2025-09-15T06:34:23Z")

</div>

from fenics import \*  
import numpy as np

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

# Parameters

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

rho = 1.0 # density  
mu = 0.01 # viscosity  
kappa = 0.01 # thermal diffusivity  
beta = 0.01 # thermal expansion  
g = Constant((0.0, -9.81)) # gravity  
T\_in = 290.0  
T\_hot = 310.0  
dt = 0.01  
t\_end = 0.5

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

# Mesh

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

nx, ny = 40, 40  
mesh = UnitSquareMesh(nx, ny)

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

# Function spaces

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

V = VectorElement(“P”, mesh.ufl\_cell(), 2) # velocity  
Q = FiniteElement(“P”, mesh.ufl\_cell(), 1) # pressure  
S = FiniteElement(“P”, mesh.ufl\_cell(), 1) # temperature  
W = FunctionSpace(mesh, MixedElement([V, Q, S]))

w = Function(W) # current solution  
w\_n = Function(W) # previous solution  
(u, p, T) = split(w)  
(u\_n, p\_n, T\_n) = split(w\_n)  
(v, q, s) = TestFunctions(W)

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

# Boundary conditions

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

inlet = “near(x[0], 0.0)”  
outlet = “near(x[0], 1.0)”  
walls = “near(x[1], 0.0) || near(x[1], 1.0)”

bcu\_inlet = DirichletBC(W.sub(0), Constant((1.0, 0.0)), inlet)  
bcu\_walls = DirichletBC(W.sub(0), Constant((0.0, 0.0)), walls)  
bcp\_outlet = DirichletBC(W.sub(1), Constant(0.0), outlet)  
bcT\_inlet = DirichletBC(W.sub(2), Constant(T\_in), inlet)  
bcT\_hotwall = DirichletBC(W.sub(2), Constant(T\_hot), “near(x[1], 0.0)”)  
bcs = [bcu\_inlet, bcu\_walls, bcp\_outlet, bcT\_inlet, bcT\_hotwall]

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

# Weak forms

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

# Time discretization

u\_mid = 0.5\*(u + u\_n)  
T\_mid = 0.5\*(T + T\_n)

# Momentum

F\_m = rho_dot((u - u\_n)/dt, v)dx   
+ rhodot(dot(u\_mid, nabla\_grad(u\_mid)), v)dx   
+ 2mu_inner(sym(grad(u\_mid)), sym(grad(v)))_dx   
- div(v)pdx   
- q_div(u)_dx   
- rho_beta\*(T\_mid - T\_in)\*dot(g, v)\*dx

# Energy

F\_T = (T - T\_n)/dt_s_dx   
+ dot(u\_mid, grad(T\_mid))_s_dx   
+ kappa\*dot(grad(T\_mid), grad(s))\*dx

F = F\_m + F\_T

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

# Stabilization (SUPG/PSPG)

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

h = CellDiameter(mesh)  
u\_mag = sqrt(dot(u\_mid, u\_mid) + 1e-10)  
tau\_m = h/(2.0_u\_mag) # momentum stabilization  
tau\_c = h_h/(2.0\*kappa) # thermal stabilization

r\_m = rho\*((u - u\_n)/dt + dot(u\_mid, nabla\_grad(u\_mid)))   
- div(2_mu_sym(grad(u\_mid))) + grad(p)   
- rho_beta_(T\_mid - T\_in)_g  
r\_T = (T - T\_n)/dt + dot(u\_mid, grad(T\_mid)) - kappa_div(grad(T\_mid))

F += tau\_m\*dot(r\_m, dot(u\_mid, grad(v)))\*dx \

- tau\_c_r\_T_dot(u\_mid, grad(s))\*dx

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

# Solve nonlinear system

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

J = derivative(F, w)

problem = NonlinearVariationalProblem(F, w, bcs, J)  
solver = NonlinearVariationalSolver(problem)

prm = solver.parameters  
prm[“nonlinear\_solver”] = “newton”  
prm[“newton\_solver”][“absolute\_tolerance”] = 1e-8  
prm[“newton\_solver”][“relative\_tolerance”] = 1e-7  
prm[“newton\_solver”][“maximum\_iterations”] = 25  
prm[“newton\_solver”][“linear\_solver”] = “mumps”

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

# Time stepping

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

t = 0.0  
vtkfile\_u = File(“velocity.pvd”)  
vtkfile\_T = File(“temperature.pvd”)

while t \< t\_end:  
t += dt  
print(f"Time = {t:.3f}")  
solver.solve()  
(u\_sol, p\_sol, T\_sol) = w.split()  
vtkfile\_u \<\< (u\_sol, t)  
vtkfile\_T \<\< (T\_sol, t)  
w\_n.assign(w)

My solver keeps getting diverging

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 15, 2025, 9:34am UTC](https://fenicsproject.discourse.group/t/2d-buoyancy-driven-flow-boussinesq-through-inlet-patch-at-the-top-wall-with-circular-obstacle-in-the-room/18210/2 "2025-09-15T09:34:09Z")

</div>

Please format you code appropriately with

````python
```python
# Add code here
```

````

and also add the output of your program as

````bash
```bash
$ Add output here
```

````

---

<div class="post-metadata">

**Author:** ![Abrar\_Nhavakar](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/abrar_nhavakar/32/8064_2.png) [@Abrar\_Nhavakar](https://fenicsproject.discourse.group/u/Abrar_Nhavakar)\
**Post date:** [September 15, 2025, 10:14am UTC](https://fenicsproject.discourse.group/t/2d-buoyancy-driven-flow-boussinesq-through-inlet-patch-at-the-top-wall-with-circular-obstacle-in-the-room/18210/3 "2025-09-15T10:14:45Z")

</div>

Python

```auto
from fenics import *
from mshr import *
import numpy as np
import os

# --------------------------
# Simulation parameters
# --------------------------
rho = 1.0
mu = 0.01
U_in = 0.01
dt = 0.001
T = 0.5
num_steps = int(T/dt)
ramp_time = 0.2

# Obstacle
obstacle_center = (0.5, 0.5)
obstacle_radius = 0.1

# Mesh refinement
base_res = 64
refine_radius = 0.15

# Arc inlet parameters
theta_start = 295 # degrees
theta_end = 360 # degrees
inlet_thickness = 0.01

# --------------------------
# Create geometry with mshr
# --------------------------
domain = Rectangle(Point(0.0, 0.0), Point(1.0, 1.0)) - Circle(Point(*obstacle_center), obstacle_radius, 64)
mesh = generate_mesh(domain, base_res)

# --------------------------
# Refine mesh near obstacle
# --------------------------
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in cells(mesh):
    x, y = cell.midpoint().x(), cell.midpoint().y()
    if (x - obstacle_center[0])**2 + (y - obstacle_center[1])**2 < refine_radius**2:
        cell_markers[cell] = True
mesh = refine(mesh, cell_markers)

# --------------------------
# Mixed function space
# --------------------------
V_elem = VectorElement("P", mesh.ufl_cell(), 2)
Q_elem = FiniteElement("P", mesh.ufl_cell(), 1)
W_elem = MixedElement([V_elem, Q_elem])
W = FunctionSpace(mesh, W_elem)
V_space = W.sub(0).collapse()
Q_space = W.sub(1).collapse()

# --------------------------
# ParaView output
# --------------------------
os.makedirs("cavity_flow", exist_ok=True)
vtkfile_u = File("cavity_flow/velocity.pvd")
vtkfile_p = File("cavity_flow/pressure.pvd")

# --------------------------
# Boundary definitions
# --------------------------
tol = 1E-14

class BottomWall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1],0.0, tol)

class LeftWall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0.0, tol)

class RightWall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],1.0, tol)

bottom_wall = BottomWall()
left_wall = LeftWall()
right_wall = RightWall()

# --------------------------
# Arc-shaped radial inlet
# --------------------------
class RadialArcInlet(UserExpression):
    def __init__ (self, center, r_in, theta_start, theta_end, U_in, **kwargs):
        super(). __init__ (**kwargs)
        self.xc, self.yc = center
        self.r_in = r_in
        self.theta_start = np.deg2rad(theta_start)
        self.theta_end = np.deg2rad(theta_end)
        self.U_in = U_in

    def eval(self, values, x):
        dx = x[0] - self.xc
        dy = x[1] - self.yc
        r = np.sqrt(dx **2 + dy** 2)
        if r == 0.0:
            values[0] = values[1] = 0.0
            return
        theta = np.arctan2(dy, dx)
        # Normalize theta to 0-2pi
        if theta < 0.0:
            theta += 2*np.pi
        if self.theta_start <= theta <= self.theta_end and abs(r - self.r_in) <= inlet_thickness:
            values[0] = self.U_in * dx / r
            values[1] = self.U_in * dy / r
        else:
            values[0] = 0.0
            values[1] = 0.0

    def value_shape(self):
        return (2,)

# Create inlet expression
inlet_expr = RadialArcInlet(obstacle_center, obstacle_radius, theta_start, theta_end, U_in, degree=2)

# --------------------------
# Boundary conditions
# --------------------------
bcu_inlet = DirichletBC(W.sub(0), inlet_expr, 'on_boundary') # Apply to entire boundary; velocity = 0 elsewhere automatically
bcu_bottom = DirichletBC(W.sub(0), Constant((0,0)), bottom_wall)
bcu_left = DirichletBC(W.sub(0).sub(0), Constant(0.0), left_wall)
bcu_right = DirichletBC(W.sub(0).sub(0), Constant(0.0), right_wall)
bcs = [bcu_inlet, bcu_bottom, bcu_left, bcu_right]

# --------------------------
# Functions
# --------------------------
w = Function(W)
w_n = Function(W)
(u, p) = split(w)
(u_n, p_n) = split(w_n)
(v, q) = TestFunctions(W)

# Initialize velocity + zero pressure
u_init = interpolate(inlet_expr, V_space)
p_init = interpolate(Constant(0.0), Q_space)
assign(w_n, [u_init, p_init])

# --------------------------
# Variational form (BDF1)
# --------------------------
F = (rho*dot((u - u_n)/dt, v)*dx
     + rho*dot(dot(u_n, nabla_grad(u)), v)*dx
     + mu*inner(grad(u), grad(v))*dx
     - div(v)*p*dx
     + q*div(u)*dx)

# --------------------------
# Time-stepping
# --------------------------
t = 0.0
for n in range(num_steps):
    t += dt
    ramp_factor = min(1.0, t/ramp_time)
    inlet_expr.U_in = U_in * ramp_factor

    # First step linearization for stability
    if n == 0:
        F_linear = (rho*dot((u - u_n)/dt, v)*dx
                    + mu*inner(grad(u), grad(v))*dx
                    - div(v)*p*dx
                    + q*div(u)*dx)
        problem = NonlinearVariationalProblem(F_linear, w, bcs, J=derivative(F_linear, w))
    else:
        problem = NonlinearVariationalProblem(F, w, bcs, J=derivative(F, w))

    solver = NonlinearVariationalSolver(problem)
    solver.parameters['newton_solver']['absolute_tolerance'] = 1E-8
    solver.parameters['newton_solver']['relative_tolerance'] = 1E-6
    solver.parameters['newton_solver']['maximum_iterations'] = 50
    solver.parameters['newton_solver']['report'] = True
    solver.parameters['newton_solver']['linear_solver'] = 'mumps'

    solver.solve()

    (u_sol, p_sol) = w.split()
    vtkfile_u << (u_sol, t)
    vtkfile_p << (p_sol, t)
    w_n.assign(w)

print("Simulation completed.")

```

Newton iteration 49: r (abs) = 5.000e-03 (tol = 1.000e-08) r (rel) = 3.768e+00 (tol = 1.000e-06)  
Newton iteration 50: r (abs) = 5.129e-03 (tol = 1.000e-08) r (rel) = 3.864e+00 (tol = 1.000e-06)  
Traceback (most recent call last):  
File “CavityFlow5.py”, line 186, in   
solver.solve()  
RuntimeError:

\*\*\* -------------------------------------------------------------------------

\*\*\* Remember to include the error message listed below and, if possible,  
\*\*\* include a _minimal_ running example to reproduce the error.

* * *

\*\*\* -------------------------------------------------------------------------  
\*\*\* Error: Unable to solve nonlinear system with NewtonSolver.  
\*\*\* Reason: Newton solver did not converge because maximum number of iterations reached.  
\*\*\* Where: This error was encountered inside NewtonSolver.cpp.  
\*\*\* Process: 0

* * *

\*\*\* DOLFIN version: 2019.1.0  
\*\*\* Git changeset: 828b3965c3ebcd128da494c596c78b17c0a303da  
\*\*\* -------------------------------------------------------------------------

```auto

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 15, 2025, 1:21pm UTC](https://fenicsproject.discourse.group/t/2d-buoyancy-driven-flow-boussinesq-through-inlet-patch-at-the-top-wall-with-circular-obstacle-in-the-room/18210/4 "2025-09-15T13:21:28Z")

</div>

A problem here is that you use `mshr`, which is not available on most systems any more.  
Have you considered switching to DOLFINx? As legacy FEniCS has not been maintained for over 5 years now.

Secondly, as your problem runs for a while, have you looked at the solution and seen if you can spot an instability building up in the proximity of the circular object?